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CHAPTER I 


INTRODUCTION 


The problem of acquisition and identification of a landmark with- 
in a given field of view is treated here from two points of view: 

1. Identification of a landmark: 

2. Estimation of its translation and rotation with respect 
to the reference frame. 

One application of this approach is to a navigation problem. One 
may have a photograph of an island and the coordinates and altitude at 
which the photograph was recorded. If at a later time a camera carrying 
vehicle flies over this island at the same altitude, then as the island 
comes into the field of view of the camera, one can, by the approach 
presented here, estimate position (translation) and orientation (rotation) 
of the craft with respect to the island. These estimates then could be 
used to command the propulsion system and navigate the vehicle. 

In addition to applications in landmark identification and acquisi- 
tion, this approach is potentially useful in problems of automatic docking 
since it permits measurement of rotation and translation of the docking 
target with respect to docking craft ( command module) . This means a 
television camera onboard the docking craft takes a picture of the docking 
target. By detection of the rotation and translation of the docking target 
with respect to the stored reference, the docking craft can position itself 
for automatic docking. 

The approach presented here requires edge enhancement 
so that the boundaries of the landmark are detected and uses the informa- 
tion contained in the boundary of the pattern by successively reading the 
coordinates of the boundary and developing a nonlinear regression analysis 
technique for simultaneous estimation of rotation and translation of the 
landmark. This method appears to be very sensitive and offers high 
resolution both in rotation and translation parameters. 



In this research only two-dimensional landmarks or patterns 
were considered. Specifically patterns in the form of ellipses and 
rectangles were first considered. The motivation for these two classes 
was to consider a class of simple shapes that can be analytically repre- 
sented, and another class that could not be analytically represented. In 
addition, different amounts of sensor noise and measurement noise were 
added to the coordinates of the boundary points to check the performance 
of this method under a variety of circumstances. 

While some alternative methods such as detection of centroid^ etc. , 
may be more useful for recognition of rectangles and ellipses, noting 
that by recognition we mean measuring rotation, translation and size of 
a pattern, it was felt necessary to consider a more general approach 
that would be applicable to more classes of two dimensional patterns 
as well as three-dimensional patterns. 

Since the method is based on tracking the boundary of a pattern, 
a review of the state of the art in boundary tracking and its application 
in pattern recognition, estimation, etc. , is provided in Chapter II under 
the title of quantization and encoding of arbitrary curves. In Chapter 
III, the fundamental nonlinear regression analysis approach is discussed. 

In Chapter IV the recognition of elliptical planar patterns is presented. 

In Chapter V the recognition of rectangular planar patterns is discussed. 
Conclusions, a summary of the results, and potential future research 
areas are discussed in Chapter VI. For use of interested readers, the 
main computer program is also documented in an appendix. 
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CHAPTER II 


QUANTIZATION AND ENCODING OF 
ARBITRARY CURVES 

This chapter is concerned with a review of approaches to boundary 
tracking and implementation of the boundary information in recognition, 
coding, estimation, etc. A more common name associated with this 
area has been contour tracing which has been used in the field of pattern 
recognition and specifically in feature extraction techniques. 

In the general problem of pattern recognition, many researchers 
feel that "contours carry a significant fraction of the information re- 
quired for recognition of image objects'^ l] . Since the recognition 
scheme developed in this research uses contour information exclusively, 
it seems appropriate to review some of the work which has been done in 
this area. 

Since most pattern recognition schemes are carried out on a 
digital computer, it is necessary to be able to represent a pattern in a 
form which may be easily manipulated by a digital computer. More 
specifically, if one is given a pictorial representation of some planar 
configuration, it is desirable to quantize and encode the boundary curve 
of this pattern into a form such that the digital computer can easily find 
such properties of the pattern as area, length of the boundary curve, 
width, height, and others to be discussed later. 

A great deal of work in quantizing and encoding arbitrary plane 
curves has been done by many researchers. [ Z -9] It is the intent of this 
chapter to review some of this research, particularly that of the chain 
representation as developed by H. Freeman. [ 10-16] 

Of the many ways in which an arbitrary planar curve ( assumed to 
be continuous) may be quantized, a particularly simple technique is 
called the grid-inter sect quantization method. In this method the curve 
is placed over a square grid, and the grid node lying closest to the 
point of intersection of the curve with a given grid line is considered to 
be a point on the quantized curve. Such a grid node is called a curve 
point. This procedure is illustrated in Figure 1, where the separation 
between adjacent grid nodes is T. 
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Fig. l--Node points for a continuous curve. 


The lines connecting adjacent curve points have length T or 
\f2 T, as seen from Figure 1. The quantized curve becomes a better 
approximation to the original curve as the grid separation, T, becomes 
small compared to the smallest instantaneous radius of curvature of 
the original curve. Freeman [ 11] points out that the grid intersect 
quantization method has an advantage over similar quantization tech- 
niques in that it comes the closest to giving equal probability to the 
occurrence of adjacent curve points which are diagonal. For an arbi- 
trary curve one would expect one half of the adjacent curve points to be 
connected by diagonal lines and one half to be connected by horizontal 
and vertical lines. 

Once the curve points have been determined, it is desirable to 
encode these points in some manner that affords economy in computer 
storage requirements and permits analytical manipulations of the pattern 
to be accomplished. One obvious encoding would be to simply store the 
coordinates of each of the curve points. However, even for a relatively 
coarse grid ( say 1024 by 1024) , each curve point would require 10 bits 
for each of its coordinates. A more economical encoding scheme takes 
advantage of the fact that since the curve which was quantized is con- 
tinuous, then successive curve points must be adjacent, as shown pic- 
torially in Figure 2. The center node is assumed to be a curve point 
and the next curve point must be one of the eight nodes shown. 

If the straight lines which join the center node with each of the 
surrounding eight nodes are assigned the same number as the correspond- 
ing outer node, then the original curve may be represented by a sequence 
of short line segments, with each line segment encoded by an integer 
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Fig. 2- -Numbering scheme for adjacent curve 
points. 



Fig. 3- -Chain representation of a continuous 
curve. 
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between 0 and 7. A line segment connecting two adjacent curve points 
is referred to as an element, and the sequence of elements which repre- 
sents the curve is called a chain. Thus, curve A in Figure 3 may be 
represented by the straight line segment curve B which is characterized 
by the chain 112221107765667. If the absolute location of curve A is 
required with respect to the x-y coordinate system, then the starting 
point, or initium, of the chain must be specified. In this case it is 
denoted by ( x 0 , y 0 ) . It is apparent that the dimension of the measure- 
ment space has been drastically reduced since now, with the exception 
of the starting point, each element ( and therefore each curve point) 
requires only 3 bits of computer storage to specify it compared to 20 bits 
which are necessary to specify the coordinates of each curve point. 
Actually, if the value of a given element is known, then the next element, 
in general, will not assume each of its eight possible values with equal 
probability. This fact can be used to further increase the coding 
efficiency by employing a chain-difference encoding scheme. [ 11] 

The chain B, which is the straight line segment representation 
of curve A in Figure 3, may be written using the chaining or "concate- 
nation" operator C defined by 


B - bj b 2 - - - b n ~ C bj 

i = l 

where b^ = 0, 1, 2, 3,4, 5, 6, or 7 

and the element b^ connects curve points i -1 and i. It is apparent that 
the number of elements in a chain will be proportional to the length of 
the curve and inversely proportional to the grid separation, T. Further- 
more, for a curve that is quantized into n curve points, the associated 
chain will have n- 1 elements if the curve is open, and n elements if the 
curve is closed. It is also readily apparent that the angle which an 
element makes with respect to the positive x-axis is simply the element 
value multiplied by 45°. 

Before considering some of the properties which chains possess, 
a few of the ambiguities in the chain representation should be pointed 
out. Consider, for instance, the chain given in Figure 3. If the abso- 
lute position of this chain is not important, then the coordinates of the 
initium, ( Xq , y 0 ) , can be disregarded and the chain is given by 
112221 107765667 as before. However, the chain may also be written as 
322123345566655. This chain represents exactly the same straight line 
segment curve, but traced in the reverse direction. Chains possessing 
this property are called inverses. Similarly, elements whose slopes 
differ by 180° are called inverse elements. Therefore, the elements 
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are inverses if 



-1 

and a^ 


-1 

a i = a^ + 4 


where the symbol + designates modulo eight addition. Thus, an inverse 
chain is obtained by finding the inverse elements of the elements of the 
original chain and then reversing the order of the inverse elements. 


From the above discussion it is apparent that any simple open 
curve (no self-intersections) has two chain representations, each being 
the inverse of the other. A simple closed curve, on the other hand, 
may be represented by any one of 2n different chains. This is due to 
the fact that there is no unique starting point for the chain; in fact, the 
chain can start at any one of the n curve points to give a total of n 
different chain representations. Some of the ambiguity in the chain 
representation of a closed curve can be eliminated if, for instance, the 
curve is always traversed in the clockwise direction, and the starting 
point is always chosen to be the curve point which is nearest to the origin. 

It is now appropriate to consider some of the properties which 
chains possess. These properties may then be incorporated into a 
variety of pattern classification schemes. It is seen that a chain is 
invariant with translation; that is to say, a chain becomes fixed in a 
coordinate system only after the coordinates of its initium have been 
specified. A chain may be rotated by k'4 5° by the modulo eight addition 
of k to each element of the chain, where k is an integer. However, the 
rotation is distortion-free only when k is an even integer, since when 
k is an odd integer the length of any element in the original chain is 
changed from T to \T 2 T or vice versa. 

The length of a chain may be directly computed by counting the 
number of even elements, n e , and the number of odd elements, n Q . 

Since an even element has length T, while an odd element has length 
\T 2 T, the length of the chain is simply 

L = n e T + n Q \T 2 T 
= ( n e + n 0 n/~ 2) T 

If n, the total number of elements of the chain, is large then the 
length of the chain may be approximated by 

L M 1 + .414 p) n T 

where p is the fraction of the adjacent curve points which are diagonal 
for the particular quantization method being used. Since p =0.41 for 
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grid-intersect quantization [ 11] , the length of a long chain is approxi- 
mately 


L % 1. 17 n T 

The height and width of a chain may also be simply computed. 
The x and y components, a x . and ay^, of each of the elements a^ are 
shown as follows: 

a i a xi a yi 

0 T 0 

ITT 

2 0 T 

3 -T T 

4 -T 0 

5 -T -T 

6 0 -T 

7 T -T 


The height is then found by subtracting the chain' s largest negative 
deviation from the x-axis from the largest positive deviation. Thus, 
the height is given by 


where 


H = (Hi) - ( Hj . 

' 17 max ' v mm 

i = 0,1,2,"*, n 


i-Z 

j = 1 


Vi 


and 


Ho E 0 


Likewise, the width of a chain is given by 

w = - < W i> 

where i = n 


mm 


w i ■ z a *i 


j =i 
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and 


W 0 E 0 


Notice that if H n = W n = 0, then the chain is closed. 

The area enclosed by a simple closed chain may also be simply 
computed. It can be shown that the area is given by [ 10, 12] 


1 1 

Area = Y a Xi ( H i_i + i a y .) 


1 = 1 


The area will be a positive number if the chain is traversed in the clock- 
wise direction, and a negative number if the chain is traversed in the 
counterclockwise direction. 

Many other properties of chains may also be determined which 
can be employed in a pattern recognition scheme. For instance, it is 
possible to determine the moments of a chain about specified axes, the 
location of a chain' s centroid, and the axes ( if any) about which a chain 
is symmetric. [12] 

Two other useful properties involve correlation functions, i. e. , 
autocorrelation and crosscorrelation. The autocorrelation function of a 

chain C a- may be defined as 
i = l 

n 

i = i 

for j = 0,+l,+2, - ‘ # ,+n 

The product a^ai + j is defined to be the cosine of the angle between 
elements a^ and a^ + j. For convenience it is assumed that the chain C 
is periodic, having a maximum period of length n. Thus 1 “ 1 

a i ~ a i + kn 

for k = 0 , +1, + 2, * * • 

The autocorrelation function is therefore defined for all j, being 
periodic with maximum period of length n. 
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The crosscorrelation function of two chains 
may be defined in two ways, 


n 

C 

i = 1 


m 

and C 
i = i 


bi 


n 

W j) = - X a i b i + j 

n t-j 


X a i + J b i 

i = i 

depending upon which chain is shifted. Here again, both chains are 
assumed to be periodic, i. e. , 

a i = a i + kn 
bi = bi 

for k = 0 , + 1, +2,*'' 

It is easy to see that the crosscorrelation function is also periodic, 
having the same period as the length of the chain which is being shifted. 
Thus , 


4>ab( J) = 4> a b( 3 + km ) 

for k = 0, + 1, I 2, * • • 

If both chains have the same length ( n = m) , then 

<f>ab(-j) = 4>ba(j) 

Since the autocorrelation function is not unique, i. e. , several 
patterns may possess the same autocorrelation function, it may only be 
used to place the unknown pattern into a class of patterns. The cross- 
correlation functions of the unknown pattern and all the patterns within 
this selected class may then be compared for recognition purposes. 
Generally the peak of each crosscorrelation function is determined, and 
recognition is based on the pattern resulting in the maximum peak. The 
crosscorrelation method has been quite effective in fitting a segment of 
a curve to a larger curve, provided that the relative scale and orientation 
are known for both curves. [ 1 5] 
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Another useful property of chains for recognition purposes is 
the so-called directionality spectrum. This consists of tabulating the 
number of elements of the chain having values 0 through 7. The direc- 
tionality spectrum is' then found by multiplying the number of odd-valued 
elements by nT 2, and drawing a bar graph of the results. A normalized 
directionality spectrum may be obtained by dividing the number of 
elements having any given value by the total number of elements 
( n e + N r 2 n 0 ) . 

A property of chains which is rotation invariant is the curvature 
property. [ 13] The curvature function of the chain & a ^ is defined 

u i + i = a i + i " a i + 8k 

where k is chosen to be - 1 , 0, or 1 so that 

I Ui + i| <4 

The sequence u^ + k is seen to be the slope change ( curvature) of the chain 

5 a, as it is traversed, 
i = i 

The number and location of the zero crossings of a smoothed 
curvature function may then be used for recognition purposes. [13] 
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CHAPTER III 


NONLINEAR REGRESSION ANALYSIS 


3. 1 Introduction 

This chapter is concerned with the basic nonlinear regression 
method of analysis that was developed for the purpose of landmark 
tracking. Since the first class’ of patterns considered are ellipses, 
in section 3. 2 representation of ellipse patterns are discussed. In 
section 3. 3 the parameter estimation problem is formulated, and its 
characteristics are delineated in sections 3. 3. 1 through 3. 3. 6. 

3. 2 Representation of the Ellipse Pattern 

The equation of an ellipse whose major and minor axes are co- 
incident with the w-z coordinate axes, as shown in Fig. 4, is given by 

g(w.z)**±+5i. - 1=0 (3.1) 

a 2 b 2 

or 

g(w,z)=e 1 w 2 + e 2 z 2 -1=0 ( 3. 2) 

where 

Cl = — , e 2 = i— ( 3. 3) 

a 2 b 2 

2a = diameter of the ellipse in the w-direction 

2b = diameter of the ellipse in the z-direction. 

If one wishes to express the equation of this ellipse with respect 
to an x-y coordinate system as shown in Fig. 5, the following transfor- 
mation holds 


w = ( x-A) cos G -f (y-B) sin 0 (3.4) 

z =s -( x-A) sin 0 + (y-B) cos 0 (3.5) 
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The equation for the ellipse in the x-y coordinate system then becomes 


F(x,y) a g(w,z) 


w = ( x-A) cos 0 + ( y-B) sin 0 
z = - (x-A) sin 0 + (y-B) cos 0 

= ej |{ x-A) cos 0 + ( y-B) sin 0] % e 2 [ -( x-A) sin 0 + ( y-B) cos 0] 2 

= [ e x cos 2 0 + e 2 sin 2 0] x 2 + [ 2( e x -e 2 ) cos 9 sin 0] xy 

+ [ e x sin 2 0 + e 2 cos 2 0] y 2 

+ [ -2e x cos 0( A cos 0 + B sin 0) -2e 2 sin 0 ( A sin 0-Bcos 0) ] x 

+ [ -2e x sin 0( A cos 0 4- B sin 0) + 2e 2 cos 0( A sin 0 - B cos 0) ] y 

4-[e 1 (Acos0+Bsin0) 2 + e 2 (A sin 0-Bcos 9) 2 ] 

- 1 (3.6) 


Eqn. ( 3. 6) contains five parameters which completely describe the 
ellipse. These parameters are e x , e 2 , A, B, and 0. One notes that 
Eqn. ( 3.6) is a nonlinear function of these parameters. 


Equation ( 3. 6) may be transformed into a linear function of a 
new set of parameters via a nonlinear transformation of the parameters. 
To this end, let the original parameters be denoted by the vector c and 
the new parameters by the vector p, i. e. , 


e l 


’pi" 

e 2 


P2 

A 

? = 

Ps 

B 


P4 

0 ' 


.P 5. 


and denote the nonlinear transformation by ip, i. e, , 

C = ijt(p) (3. 8) 


In order to derive the nonlinear transformation of Eqn. ( 3. 8) one may 
rewrite Eqn. ( 3. 6) as 

F( x, y) = p j x 2 + p 2 xy + p 3 y 2 + p 4 x + p 5 y + p 6 -1=0 ( 3. 9) 


where 
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X 


Fig. 5 --Rotated and translated ellipse. 



Fig. 6 --Relation of 0 and p . 
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p 1 

= ej cos 

2 e 

+ e 2 

sin 2 

e 




( 3. 10) 

P 2 

= 2( ei - 

e 2 ) 

cos 

0 sin 

e 




(3. 11) 

P 3 

= e x sin 4 

2 e 

+ e 2 

2 

cos 

e 




(3. 12) 

P 4 

= - 2A( e 

i! cos 2 0 + e 2 

sin 2 0 ) -2B( e 

x -e 2 ) cos 

0 sin 0 

(3. 13) 

P 5 

= -2B( e 1 

sin 2 0 + e 2 cos 

2 0) -2A( e 1 

-e 2 ) cos 0 

sin 0 

(3. 14) 

P6 

= A 2 ( ei 

cos 

2 e + 

e 2 sin 2 

0) + B 2 ( ei 

sin 2 0 + e 

2 cos 2 

6) 


+ 2AB( ej -e 2 ) cos 0 sin 0 (3.15) 

Eqn. ( 3. 10) -( 3. 15) may be manipulated to obtain the values of the 
parameters e x , e 2 , A, B, and 0 in terms of the parameters p j , p 2 , p 5 , 
p 4 , p 5 , and p 6 . From Eqn. ( 3. 10) , ( 3. 1 1) , and (3. 13) 

p 4 = -2Ap!-Bp 2 (3. 16) 

while Eqn. (3.11), (3.12), and (3. 14) yield 

P 5 = " A p 2 “ 2 B p j (3. 17) 

Solving Eqn. ( 3. l6) and ( 3. 17) simultaneously for A and B gives 


A = ~ ^ P 5 P 4 + P 2 P 5 

. 2 

4 P 1 P 3 ‘ P 2 

"2 Pi P 5 + P 2 P 4 

B = — 

. 2 

4 P 1 P 3 - P 2 

From Eqn. ( 3. 10) and (3. 12) one obtains 

Pi + P 3 = e i + e 2 

and pi - p 3 =(e 1 -e 2 )(cos 2 0 - sin 2 0 ) 

= ( ej -e 2 ) cos 2 0 


(3. 18) 


(3. 19) 


( 3. 20) 


(3.21) 
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while Eqn. ( 3. 11) yields 


P 2 = (e x -e 2 ) sin Z 0 ( 3. 22) 

Eqn. ( 3. 21) and ( 3. 22) then give 

tan 2 0 = — — 

Pi-Ps (3.23) 

or 

0 = | tan" 1 .. P2 ( 3. 24) 

P l “P s 


If a right triangle having s ides of length p 2 and p x - pj is formed, then the 
hypoteneuse has a length ^ p 2 2 + ( P l “P 3 ) 2 as shown in Fig. 6. From 
Fig. 6 it is apparent that 


sin 2 0 = ■ ■■■ 

n/ p 2 2 + ( Pi - P 3 ) 2 


( 3. 25) 


Then from Eqn. ( 3. 22) and ( 3. 25) one obtains 


^P2 2 + (pl"Ps) 2 “ e l “ e 2 
and then Eqn. ( 3. 20) and ( 3. 26) yield 


(3. 26 ) 


e i = 2 ( P 1 + P 3 + ^ P 2 2 + ( P 1 - P 3) 2 ( 3 - 27 ) 

e 2 = i ( P 1 + P 5 " ^ P 2 2 + ( P 1 "P 5 ) 2 ) ( 3 . 28 ) 

Since only five independent parameters are required to fully 
specify an ellipse, it is reasonable to expect that Eqn. ( 3. 9) may be 
simplified. Dividing Eqn. ( 3. 9) by ( p$ -1) gives 

F( x, y) = p ! x 2 + P 2 xy + p 3 y 2 + p 4 x + p 5 y + 1 = 0 ( 3. 29) 

where 

p t = 9 1 for i = 1, Z, 3,4, 5 (3.30) 

P 6 -1 
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r~ 

From Eqn. ( 3. 18) , ( 3. 19) , and ( 3. 24) one observes that A, B, and 0 
are ratios of the M p M parameters where both numerators and denomina- 
tors are of the same order. Thus, the denominator of Eqn. ( 3. 30) will 
cancel, making Eqn. ( 3. 18) , ( 3. 19) , and ( 3. 24) the same function of 
the M p M parameters as of the "p 11 parameters. Therefore 

A = .-gPiP. 4 + -P 2 P 5 

4piPs-P2 2 (3.31) 

B = .-gP.lP5 +P2P4 

4piPs-Pz 2 (3.32) 

0 = \ tan -1 -2l 

Pi "Ps 

( 3. 33) 

This is not the case for Eqn. ( 3. 27) and ( 3. 28) , however. They become 
e i = (p6-!)[ i(pi + Ps + ^ P 2 * +( Pi -Ps) 2 ) ] ( 3. 34) 

e 2 = (p6-l)ti(Pi + Ps - P2 2 + (Pi -Ps) 2 ) 1 (3.35) 

From Eqn. ( 3. 10) , ( 3. 1 2) , and ( 3. 15) 

P6 = A p i + B p j + ABp j (3, 36) 

Eqn. ( 3. 30) and (3 . 36) then give 

P 6 ~ 1 = (A 2 pi + B 2 Ps + AB p 2 ) -1 

= [( P6-1)(A 2 Pl + B 2 p s +AB p 2 )]-l (3.37) 

or 


A 2 Pj + B 2 p s + ABp 2 -1 (3.38) 

Making use of Eqn. ( 3. 31) and ( 3. 32) further reduces Eqn. ( 3. 38) to 

P6-l = I 

PlPs 2 ± Ps P4 2 -P2P4P5 - 1 

4 Pi P 3 -P 2 2 
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= 4 PtP, -p 2 2 

PiPs 2 +PsP4 2 + P 2 2 -P 2 P 4 P 5 -4piPs ^ 3 39 

By substituting Eqn. ( 3. 39) into Eqn. ( 3. 34) and ( 3. 35) the parameters 
e l and e 2 become functions only of the M p M parameters. 


4 Pi Ps -pz . 2 [_i(pi 

PlP5 2 + P3P4 2 + P2 2 "P2P4P5 “4pi Pj 


+ Ps + P 2 2 + ( Pi -Ps ) 2 ‘) ] 

(3.40) 


e 2 


4 P i Ps ~P2 . Z [i( Pl 

P 1 P 5 2 + PsP 4 2 + P 2 2 -P 2 P 4 P 5 - 4 pi Ps 


+ Ps P 2 2 + ( Pi -Ps) 2 ' ) ] 

(3.41) 


If the number found from computing Eqn. ( 3. 39) is negative then the 
expressions for e r and e 2 as given by Eqn. (3.40) and 3.41), respec- 
tively, should be interchanged as seen by referring to Eqn. ( 3. 21) , ( 3. 22) , 
and ( 3. 24 ) . 


The derivation of the transformation 


C = $ p) = 


4*i (p) 
4z (p) 
4s ( p) 
44 (p) 
4s (p)_ 


(3.4 2) 


has now been completed, with the components ij^ , 1 fj z , ijj 3 , 414 , and 1^5 given 
by Eqn. ( 3.40) , ( 3.4 1) , ( 3. 31) , ( 3. 32) , and ( 3. 33) , respectively. 


It was shown that an ellipse may be expressed as a linear func- 
tion of a set of five parameters or as a nonlinear function of another set 
of five parameters as given by Eqn. ( 3. 29) and ( 3.6) , respectively. 

The two sets of parameters are related by the transformation given by 
Eqn. ( 3. 8). It is now appropriate to investigate methods whereby the 
unknown parameter vector for an ellipse may be estimated after points 
on the ellipse have been measured. 
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3. 3 Parameter Estimation Problem 


3. 3. 1 The Error Formulation 

In order to estimate the five parameters that represent the size 
and position of the ellipse one may write Eqn. ( 3. 29) or ( 3. 6) as 

F ( x i» Yi5 To) = 0 ( 3. 43) 

where (x^,y^) is any point on the ellipse and 

ir 

£o = £s 

£4 

y (3.44) 

is the parameter vector which characterizes the ellipse. 

If. Eqn. ( 3. 29) is used for describing the ellipse, then 

Pi' 

* P2 

£o = Ps 

P4 

. PS . 

If one measures any point on this ellipse and computes F using 
some other parameter vector, C, the value for F will not be zero, but 
rather it will be equal to an error, € . That is, 

F ( x i> Yi : C) =€ (3. 45) 

Likewise, if one makes an error in the measurement of a point on this 
ellipse, then the computed value for F using the true parameter vector , 

C 0 , is again non-zero, representing an error, £ . That is, 

F( x i> Yi ; T) = E ( 3.4 6) 

'v a, 

where (x^,y^) is now a noisy measurement point. 

Now as the boundary of the landmark or the pattern is traced, a 
sequence of coordinates x^, y^ become available. Given Eqn. (3.4 5) 
and the coordinates of the boundary ( x^, y^) , i = 1, 2, " * , one can 
estimate a vector C e of the true parameter vector C 0 . 
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The estimation will be based on minimizing an appropriate func- 
tion of the error in Eqn. ( 3.4 5)* The simplest of these functions appears 
to be the sum- squared of the error. If N points on the boundary are 
available, the sum-squared error is given by 

<f>(x, y; f) = ^ [F(xi,y i ; X) - F(x i ,y i ; |^) ] 2 

i = l 
N 

\ _ »-o -*.2 

= t F ( x i’ > t] 

i = l 

For convenience, the tilde on xi and y^ will be eliminated from now on. 

It will be understood that ( x{, y^) represent noisy measurement points. 

Let 


(3.47) 


Min <f> (x, y" ; t) = <K x, yX e ) (3.48) 

r 


where x, y are N-dimensional vectors consisting of the N points which 
were measured on the ellipse or landmark boundary. 


By defining an N x 1 error vector, e 


e = F(x, y ; £) 


F(x i ,y i ; £) 

-F( X N’ ^ 1 


the sum-squared error may be conveniently written as 


<Kx, y ; X ) 


-T - 

e e 


(3.49) 
( 3. 50) 


where T denotes transposition. 


It should be noted that the criterion function, cj), corresponding 
to the sum-squared error is dependent upon whether F given in Eqn. (3. 49) 
corresponds to Eqn. ( 3.6) or to Eqn. ( 3. 29) . The resulting criterion 
functions are not identical. This point is discussed in greater detail in 
Appendix II. 
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3. 3.2 Error Minimization by Linear 
Regression Analysis 

When the ellipse is given by Eqn. ( 3. 29) then it is a linear func- 
tion of the parameter vector, p, and the minimization of the sum- squared 
error reduces to a simple result. The error may be written as 

e = F( x, -y ; p) = M p + T ( 3. 51) 

where M is the N x 5 matrix 


2 

X i 

XiYi 

Yi 2 

x i Yi" 

K 

I 

X N^N 

2 

YN 

X N Y N _ 


( 3. 52) 


and 1 is an N x 1 vector containing all 1' s. Eqn. ( 3. 52) indicates that 
the elements of M are simply functions of the measured points on the 
ellipse. 


The sum-squared error, becomes 

^ ^ _*T _ 

<Kx, y ; p) = e e 

= (M p + T) T (M p +f) (3.53) 

which is a positive definite quadratic form in the coordinates of the trial 
parameter vector, p*. Therefore, one merely needs to compute all of 
the partial derivatives of (j> with respect to the components of p and equate 
them to zero to find the unique minimum value of <(>. Expanding Eqn. 

( 3. 53) gives 

cf> = ( M p) T ( M p) + (Mp) T l+l T Mp + t T l 

= p^ M p + 2 p' T + N (3. 54) 

and therefore 

|4r = 74> = 2 M T M p“ + 2 M T T . ( 3. 55) 
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and 



9j> 

T 

= 0 = 2 M Mp e + 2 M T T 



9? 


(3.56) 



P Pe 


or 

m t 

M p^ = -M T T 

(3. 57) 


T 

The matrix M M may be inverted, assuming M is of full rank, to give 

P e = - [M T M] _1 M T T ( 3. 58) 


This estimate is then the M least squares estimate 11 of the true parameter 
vector, , and shall be referred to as the one step minimization method. 

It should be noted that the simple expression for p e as given in 
Eqn. ( 3. 58) would not have resulted had the criterion function been 
something other than quadratic in the parameter vector components, 
since the differentiation would have yielded a nonlinear relation in the 
parameter vector components. 

If the parameter vector, cj> , is to be estimated directly from 
Eqn. ( 3.6) some other technique than Eqn. ( 3. 58) must be employed 
since the parameters enter Eqn. ( 3. 6) in a nonlinear manner, making 
the criterion function, <\>, no longer quadratic in the parameters. A 
complete automatic computer algorithm to estimate parameter vectors 
for this sort of problem has been developed by R. B. McGhee [ 17] and shall 
be utilized here. This is an iterative minimization scheme rather than a 
one step minimization scheme such as was associated with the linear 
regression analysis. The essence of this computer algorithm is dis- 
cussed below. 

3. 3. 3 Gauss -Newton Iteration 

If the nonlinear response vector, F, has its Taylor series ex- 
pansion truncated after the linear term 


or 


F(c 1 + A "c) 


F(^) + 9_F_ 

a c 


A c 


C a Ci 


A 

F( ci + Ac) = F( Cj ) + Z A c 


(3. 59) 
(3.60) 
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where 



9 F( X] , yi 



3 Cj 


9 F( x N , y N ; c) 
3 Cj 


9 F( x, , y, ; c) 

3 C 5 


9 F(x n , y N ; c) 

3 C 5 

C = Ci 

(3.61) 


then the criterion function, cj>, associated with this response function is 

A A rji A 

4>( x, y ; A c) c\ ) = F F (3.62) 


A 

where F is defined in Eqn. ( 3. 60) and 

a T 

<|)( A c) = ( cT + Z A c) ( e + Z A c) (3. 63) 

which is a quadratic form in A?. Expanding Eqn. (3. 63) gives 

<j)( A c) = e^^e + e^ZA c + A c^Z^e + A cT^Z^Z Ac (3.64) 

If Eqn. (3. 64) i s differentiated with respect to Ac and the result equated 
to zero, the minimizing value of A c becomes 

T + T 

= 0= 2 Z e + 2 Z ZAcj 


3 

3 Ac 


A cT = A c x 


(3.65) 


or 

T -l T 

A ci = - ( Z Z) Z e (3. 66) 

T 

assuming that Z Z is nonsingular. The matrix 

S = Z T Z (3. 67) 

is referred to as the regression matrix due to the similarity of this 
method to linear regression analysis. The normal equation for iteration, 
Eqn. (3. 65) > is linear in Ac only because function F was linearized and 
a quadratic criterion function was chosen. 
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Since 


f<t>( Cl ) a |4- 

9 C 


= !_? T ? 
9 c 


= Z(||) T e = 2 Z 
c = c, 18 c/ c = c, 


C = Cl 


Eqn. (3. 66) becomes 

Aci = - i s ’ 1 F<k<?i) = p; 

The new value for the parameter vector, c, then is 

c 2 = cTi + A ci 


c = c x 

( 3 . 68 ) 

(3. 69) 

(3. 70) 


upon which a new iteration may then be initiated. This procedure is 
referred to as the "Gaus s - Newton” iteration method [ 17] . 

As mentioned previously, Eqn. (3. 69) is based on the assump- 
tion that linearizing function F, Eqn. (3. 59) > is valid. Since, in fact, 
this may be completely invalid, it is quite possible that the sequence of 
parameter vector estimates, given by Eqn. (3. 70) , will not converge to 
c 0 . The necessary and sufficient conditions for the convergence of the 
Gauss-Newton' procedure may be derived; however, the test is generally 
complicated enough that in practice one simply computes <|)( cT- + f^) at 
each step to see if an improvement results. It can be shown that the 
Gauss-Newton iteration always converges when binary scale factor ad- 
justment is used [ 17] . When this technique is used, ?i + i is found 
from 

Ci + ! = Ci + 2" k Pi (3. 71) 

where k is the first non-negative integer which reduces cj>. However, 
experimental results show that the rate of convergence can be quite slow. 
For this reason, the "modified 11 Gauss -Newton procedure will not be used. 
The Gauss -Newton iteration enjoys its greatest success as a terminal 
iterative technique, where the current parameter vector is "close" to its 
minimizing value . 

3. 3.4 Newton-Baphson Iteration 

When the Gauss -Newton iteration fails to give a reduced value 
for the criterion function, <J), then direct gradient techniques may be 
appropriate. This eliminates the necessity of inverting the matrix, S, 
and also makes it possible to handle parameter range constraints in a 
straightforward manner. 
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The gradient technique to be employed is the method of steepest 
descent, in which case the parameter change vector, Ac, is directly 
proportional to the negative gradient of the criterion function, <|>. 

Ac t = - k f<Kci) (3.72) 

where 

k > 0 (3.73) 

Thus, Ac{ is in the direction of the greatest rate of decrease of <{>. The 
next parameter vector estimate then becomes 

c i + i = c i * ^ c i ( 3. 74) 

which may then be used to perform another iteration. 

Before Eqn. ( 3. 72) can be utilized, it is necessary to choose 
some value for the scale factor, k. The n Newton-R aphson 11 method may 
be used to obtain a value for k. Essentially it is based on taking the 
linear portion of the Taylor series expansion of the criterion function, cj>, 
and extrapolating this to zero. More precisely, suppose that (j> is a 
sufficiently smooth function of ? such that it may be represented locally 
by the Taylor series 

<KcT+ Ac) a <f<c) + F<j> T Ac + 0( A? 2 ) (3.75) 


^2 

where 0(A c ) represents all the terms in the series which are quadratic 
or higher order in Ac. Then for small Ac, 0( Ac? 2 ) may be ignored, and 



^ rp ^ 

cj>( c + A c ) = cj)( cT) Hh V <\> Ac 

(3.76) 

or 

<(>(?+ Ac) = cj>( *c) - k |^^| 2 

(3.77) 

using Eqn. ( 3. 72) . 

0 

Extrapolating <|> to zero, for which k = k , 

gives 


0 = <K Si) -k° |p<K Si) | 2 

( 3.78) 

or 

o 4>( q ) 

kl I^Ci)! 2 

(3.79) 
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The corresponding parameter change vector, Aq, at each state of itera- 
tion is then, from Eqn. (3. 72) , 


A 


-J) 

c i 


tfci) m ci) 
P # cj | 2 


( 3. 80) 


It may well be the case that iteration based on Eqn. ( 3. 80) will 
give a value for c* + A c) which is larger than cj>( cT) . This simply means 
that Ac*is so large that linear extrapolation of <\> to zero is, in fact, 
invalid. Eqn. (3. 75) guarantees that for some 0 < k < k° the 
criterion function, cj>, will be reduced, however. Thus it is desired to 
find some k = such that 

# 

Min<Kc'i-k i ^(j)( ci))= <j)( Ci'ki 7<Kci)) 

k > 0 ( 3. 81) 


The next parameter vector estimate is then 
Ci + 1 = Ci-ki r 4>( Ci) 


( 3. 82) 


This iteration scheme is called the "optimum gradient method 11 . 


It is, of course, not feasible to search over all values of k on a 
computer, but it is quite feasible to perform a binary search over the 
range 

o < k < f = k° (3.83) 

I Fc|> | 

which may be considered a "suboptimum gradient method". Assuming 
that <\> is continuous and ff<J> ^ (T, Eqn. ( 3. 7 5) guarantees that there exists 
an n such that 


<j>( ?i- J-n k? f§) < 4< cT.) 
2 


(3.84) 


and therefore a binary search procedure always produces a convergent 
sequence of values for cj>. A simple algorithm may be constructed to 
find the minimizing value for n as follows. First of all, compute Ac^ 
from Eqn. (3. 80) . Then, for n = 0, 1, ' ' ' , evaluate 


, n ,/ -*• 

<t>i = 4>( ci + 


l 

2 n 



( 3.85) 
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Once a value for n is reached, say n = m, such that 


,m-i . ,m ^ . m + i , 0 x 

9i > 9i < 4>i ( 3. 86) 

and 

<*i ( 3 * s? ) 

then take for the new value of c 

c i + i = ^ + _L Ac? 

2 (3.88) 

If Eqn. ( 3. 87) is not satisfied, continue increasing n until Eqn. ( 3. 86) is 
again satisfied and then check Eqn. ( 3. 87) once again. Continue this 
until both Eqn. ( 3. 86) and ( 3 .87) are satisfied. Eqn. ( 3. 88) is then 
the new value for c. 

A further refinement may be incorporated into this algorithm by 
fitting a quadratic function to the points cf)^ and + 1 . Letting 

= 4b , = <f>i and 1 = 4 * 2 * it is straightforward to show 

that the minimum of this fitted quadratic function occurs at 


where 


If 



4*2 - 5<t>i + 4<j?o 
4*2 -34*i + 2 4 *d 



k° 

i 


4>( c i-kq ua d J^4*) < 4*( c i - ^ 1 i ^4*) 


then the parameter change vector is 

A ct = - k qua d V<\>( cO 

Otherwise, 

A ci = - kj ?.) 

3. 3. 5 Pange Constraints 


( 3. 89) 


( 3. 90) 


( 3.91) 


(3.92) 


It is necessary to place constraints on the parameters since we 
assume the landmark or ellipse to be in the field of view. These are 
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range constraints, in which each component of the parameter vector is 
independently restricted to lie within some specified interval on the 
number scale. Thus, each component, c^, must satisfy 

a i £ ci £ b^ ( 3. 93) 

where a^ and b^ are the lower and upper limits of the allowed range, 
respectively. This then means that the parameter vector, c, which 
minimizes the criterion function, cj>, must be in a hypercube in parameter 
space. 


Considering the problem at hand, one notes that in order for Eqn. 

( 3. 2) to represent an ellipse it is necessary that and e 2 ( or c l and 
c 2 , respectively) be positive. Likewise, the finite field of view of the 
optical equipment places constraints on e 2 and e 2 as well as the trans- 
lations A and B ( or c 3 and c 4 , respectively). Since an ellipse is 
symmetric about its two axes, the rotation angle, 0 = c 5 , may be con- 
strained to lie in the first quadrant. 

Since the gradient-descent method discussed earlier is valid only 
on the interior of the 5-dimensional contstraint region, R, it is necessary 
to use a different strategy when a constraint boundary is encountered 
during a gradient-descent. The method to be used is called the gradient- 
projection method. If a constraint boundary should be encountered, this 
method projects the gradient onto the constraint surface and then travels 
in the negative direction of the projected gradient until a minimum for 
<j> is found. The actual minimum for <j> may either be located on the 
interior of R or on a constraint boundary of R . In the latter case the 
projection of the gradient will have all of its components equal to zero 
at the final iteration. 

The mechanization of the gradient-projection method may be 
performed in three steps. 

1. Check each component of the current parameter vector 
estimate, c , to see whether it is within the allowed range 
or if it lies at the lower or upper end of the range. 

2. If any component lies on either extreme of its range, and 
if the negative of the corresponding gradient component 
points out of the constraint region, then set this component 
of the gradient equal to zero. Leave all other components 
of the gradient at their true value. 

— J 

3. The resulting vector, V <\>p > is the desired projected 
gradient. 
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In order to find the optimum step size, it is necessary to find the maxi- 
mum scale factor which can be applied to - without violating a 

range constraint. To this end, suppose that is positive. This 

means that cj can be reduced in value without J violating a range 

constraint. Let k!j be the largest scale factor that can be applied to the 
negative j** 1 gradient component without violating the j** 1 range constraint. 
Then k<J satisfies 


9 d> t-v 

c ; -k? P = ai 


' J J a cj 


(3.94) 


which gives 


c i -a i 
vO _ J J 

j » 4»p 


(3.95) 


Likewise, for the negative components of 


it follows that 


c -i “b-j 


9 <t>p 
9 c-i 


(3.96) 


The maximum scale factor, k 0 , is found from 

k 0 = Min {k9 } ( 3. 9?) 

j J 

where the k<j are defined only for the non-zero components of F<t>p- 

The maximum step size for the parameter change vector now 
becomes 


4c i= - k 0 V4p (3.98) 

This value may then be used as the maximum step size for the binary 
search procedure discussed earlier. 
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3. 3.6 Global Optima 


Both the Gauss-Newton method and the Newton- Raphs on method 
are suited for determining local minima since they make use only of 
local information. If more than one minimum is contained within the 
constraint region, R, then it is desirable to find the smallest of all of 
these minima. Such a minimum is called a global minimum. Of course, 
the only way to find the global minimum with certainty is to exhaustively 
search the entire constrained parameter space. Since this is not feasible 
or practical on a computer, one must choose some method whereby a 
given confidence level is attained that the minimizing parameter vector 
obtained is associated with a value of cf) which is smaller than some 
specified per cent of the points in R. Such a method is that of uniform 
random searching in which parameter vectors are chosen at random 
(with a uniform distribution for each component) and their correspond- 
ing criterion functions are evaluated. The parameter vector, c, being 
associated with the smallest value for cf) is then used to initiate a local 
minimization [ 17] . 
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CHAPTER IY 


RECOGNITION OF ELLIPTICAL PLANAR PATTERNS 


4. 1 Introduction 

This chapter is concerned with the recognition of elliptical planar 
patterns. By the term ’'recognition 11 it is meant that the two minimization 
techniques which are discussed in Chapter III are employed to estimate 
the five parameters associated with an ellipse having arbitrary size and 
shape, as well as arbitrary position (translation and rotation) in the 
planar field of view. 

Section 4. 2 discusses the statement of the problem and the general 
approach which is to be pursued, while section 4. 3 discusses the vari- 
ous parameters which are associated with the implementation of the two 
minimization schemes. 

The results which were obtained from the two minimization 
schemes are discussed in section 4. 4, and a summary of the advantages 
and disadvantages of the two methods is contained in section 4. 5. 

4. 2 Statement of Problem 

The parameter estimation schemes were first applied to the 
recognition of elliptical patterns. Elliptical patterns were selected 
first because they are a somewhat complex pattern and yet their boundary 
may be represented analytically. Furthermore, most of the ground work 
for the estimation of the parameters of an ellipse has been laid in 
Chapter III. 

As was discussed in Chapter III, an ellipse, located in a plane, 
may be fully characterized by five parameters. Two parameters, e x 
and e 2 , are neceasary to specify the size and shape of an ellipse, while 
three parameters are necessary to specify its position and orientation 
in the plane. Figure 7 shows a typical ellipse which has been translated 
and rotated with respect to the reference x, y-coordinate system. 
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Fig, 7--Parameters of an ellipse in the x, y-reference 
frame. 

With respect to the w> z-coordinate system, this ellipse may be expressed 
analytically by Eqn. (4.1) 


e x w 2 + e 2 z 2 = 1 


where 


e i 



and 



(4. 1) 


(4. 2) 


The parameters r w and r z are respectively called the w-axis radius and 
the z-axis radius of the ellipse. The parameters which are actually 
estimated are e x and e 2 , which are related to r w and r z by Eqn. (4. 2). 

The x and y-translation parameters are denoted by A and B, respectively, 
and the rotation parameter is denoted by 0. These parameters are all 
shown in Figure 7. 

The parameter vector, "c , which characterizes an ellipse is given 
by Eqn. ( 4. 3) . 


’Cl" 


‘ e l 

C 2 


e 2 

C 3 

= 

A 

C* 


B 

• C 5- 


.9. 


It is, then, the intent of this chapter to determine the feasibility 
of recognizing an elliptical planar pattern by estimating its associated 
parameter vector c?, and to determine whether the one step minimization 

method or the iterative minimization scheme does the better job of per- 
forming this parameter estimation task. 
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4. 3 Implementation of the Parameter 
Estimation Schemes 


To simulate an ellipse in the field of view, points which lie on 
the boundary of an ellipse are artificially generated by the subroutine 
denoted by DATA. The logic by which this subroutine selects the data 
points is discussed in Appendix I. After the subroutine DATA is pro- 
vided with parameter vector cT 0 , it generates data points which lie on 
the boundary of an ellipse which is characterized by cf 0 . This parameter 
vector was arbitrarily chosen to be 


1 . 00 
0. 25 
1.00 
- 2 . 00 
0. 50 


(4. 4) 


This corresponds to an ellipse which has a w-radius and a z-radius equal 
to 1. 0 and 2. 0, respectively. In addition, the ellipse has been translated 
one unit in the x-direction and two units in the negative y-direction, and 
rotated 0. 5 radians. 


Thus, in the absence of measurement noise, one would expect 
the estimate for <? to be exactly c Q . 

While the parameter vector <T 0 was held fixed, two other parameters 
were varied to determine their effect on the accuracy of the parameter 
estimation schemes. 


One of these variable parameters was the number of data points 
which were used to represent the boundary of the ellipse. Ten data 
points were chosen for a sparse distribution of points on the boundary, 
while one hundred data points were chosen for a dense distribution of points 
on the boundary. Intermediate values for the number of data points were 
chosen as 20 and 50. 

The other variable parameter was the amount of noise which was 
added to the data points to simulate the effect of measurement noise or 
other errors. The noise samples, which are generated on the digital 
computer, have a gaussian distribution. The mean and standard deviation 
of these noise samples may be independently specified. In all cases the 
mean was chosen to be zero, while the standard deviation was either 
0. 0, 0. 1, 0. 2, 0. 3, 0.4, or 0. 5. The maximum value for the standard 
deviation of the noise sample, 0. 5, was one -half of the z-radius of the 
noiseless ellipse. Noise samples having a standard deviation larger than 0. 5 
result in the data points having such a large scatter that they no longer even 
remotely resemble the boundary of an ellipse. In fact, physical systems 
which correspond to the higher values of standard deviation ( 0. 3-0. 5) 
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would have limited practical utility, but it is of interest to investigate 
the reliability of the parameter estimation schemes for high noise levels, 
and to develop bounds on the performance of such systems. 

The pattern which is to be recognized is required to lie within 
some bounds since in a physical situation the optical system would have 
a finite field of view. The field of view was arbitrarily chosen to be a 
square measuring eight units on a side. The boundary of the ellipse 
which is characterized by the parameter vector c Q given by Eqn. ( 4. 4) 
lies entirely within this field of view. 

A remark should be made at this point. If the noise which is added 
to the data points has a large standard deviation, it is possible that some 
of the resulting noisy data points will fall outside of the field of view. 

When this situation arises, those noisy data points which fall outside of 
the field of view are still regarded as valid data points in the simulation. 
Physically, in an actual landmark tracking or automatic docking situation, 
noise may be classified into two general categories. The first category 
consists of noise associated with measurement errors. These include 
grid quantization errors, detector or sensor errors, and transmission 
errors. In any of these cases the coordinates of a data point (which is in 
the optical field of view) will be in error, and if the true data point is near 
the boundary of the field of view then it is possible that the noisy, or 
measured, data point will have coordinates which lie outside of the field 
of view. This situation is contrasted to the second-category into which 
noise may be classified, which may be termed "masking’ 1 noise for lack 
of a better name. This kind of noise corresponds to a case in which the 
field of view is partially covered with clouds or to a case in which the 
optical system is very badly out of focus. The sensor will be unable to 
detect the data points which are masked, or obscured, due to either of 
these situations, and therefore these data points are in essence, discarded. 
This masking noise has the effect, therefore, of shrinking the field of 
view. 


Thus, it can be seen that the noise which is being simulated cor- 
responds to measurement noise rather than "masking” noise. 

The numerical values utilized in simulation experiments for the 
range constraints for the five ellipse parameters are as follows: 


0. 0625 < e x < 16. 0 
0. 0625 < e 2 < 16 . 0 
-4. 0 < A < 4. 0 

-4. 0 < B < 4. 0 

0 < 0 < 1 . 57 


(4. 5) 
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These range constraints permit the fitted ellipse to have either of its 
radii range in size from 1/4 unit to 4 units ( i. e. , the maximum diameter 
is constrained to be no larger than the dimensions of the field of view) . 

In addition, the center of the fitted ellipse is permitted to lie anywhere 
within the field of view, while the rotation angle is constrained to lie in 
the first quadrant due to the symmetry of an ellipse. 


The estimation process has to be started with an arbitrary initial 
parameter vector, c e - The initial guess for the parameter vector was 




' 0. 5" 
1. 3 
0. 7 
-2. 5 
0 . 8 


( 4. 6) 


which corresponds to an ellipse having a w-axis radius and a z-axis 
radius equal to 1. 414 and 0. 877, respectively. The parameter vector cT e 
was chosen such that its components were "close 11 in value to the com- 
ponents of c 0 and yet not so close as to make the estimation problem 
trivial. 


After performing a local minimization using <? e as the initial 
estimate for ? Q , four more local minimizations are executed with the 
initial estimate in each case being found by the RANSEE ( random search) 
subroutine. [ 17] Thus, a total of five trial local minimizations are 
carried out. The number of random searches for each trial local minimi- 
zation was set equal to 100. 

4. 4 Results 

Both the one step minimization method and the iterative minimi- 
zation scheme which are outlined in Chapter III were employed to esti- 
mate the parameters of the given ellipse. The results which were obtained 
by using these two schemes are shown in Tables 1 and Z, respectively. 

It should be pointed out that the data points are exactly the same for both 
minimization schemes, permitting a meaningful comparison to be made. 
The results which are tabulated in Table Z are also shown pictorially in 
Figures 8, 9, 10, 11, 1Z, and 13. In these figures the ellipse having 
a solid line boundary corresponds to the parameter vector <? 0 . The "+" 
symbols correspond to the noisy data points arising from the solid line 
boundary. The ellipse which is fitted to these noisy data points, and 
characterized by c e , is represented by the dashed line boundary. The 
x and y-radii correspond to the w-axis and z-axis radii, respectively. 
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A comparison of Tables 1 and 2 shows that for noise levels below 
<r = 0. 4 the two minimization schemes produced results which were quite 
similar. Referring to these tables or to Figure 8 one notes that for 
noiseless data points the parameter vector is estimated precisely, that 
is, c e = c Q . This is a criterion which any good recognition scheme 
should fulfill, of course. 

Figure 9 shows the results which were obtained for c r = 0. 1. 

Special note should be made concerning the accuracy with which the 
rotation component of the parameter vector was estimated. It is seen 
that the largest error is less than 3 degrees, while for three of the four 
cases this error is considerably less than one degree. 

A brief comment concerning the expression for the error should 
be made at this point. In most instances it is more convenient and 
meaningful to express an error in percentage rather than absolute terms. 
Such is the case here. Since an ellipse is symmetric about both its 
vertical and horizontal axes, its angular position is unique only in the 
first quadrant, i. e. , 0 to 90 degrees. The percentage error may then 
be defined as the ratio of the absolute error to 90 degrees. With the 
percentage error so defined, one can see that for cr = 0. 1 the maximum 
error is approximately three percent for the rotation parameter, which 
is quite good considering the fact that the reference rotation angle is not 
constrained to be small. 

The results for cr = 0.2 are shown in Figure 10. Here again one 
notes that the error in estimating the rotation parameter is quite good. 

In fact, ignoring the 10 data point case, the maximum error is still less 
than three percent. In the 10 data point case the error is approximately 
seven percent, which is still reasonable considering the scarcity of data 
points and the noise level. Another observation which can be made from 
both Figures 9 and 10 is that the estimates for the parameter vector become 
better as more data points are used, a situation which intuitively seems 
reasonable. 

When the noise level reaches cr = 0. 3, as shown in Figure 11, 
the fitted ellipses begin to differ from the reference ellipses to a larger, 
and perhaps unacceptable, extent. The scatter of the data points is such 
that an accurate fit cannot be realized by either of the parameter esti- 
mation schemes. However, it should be pointed out that the estimate for 
the rotation parameter is still respectable, except for the 10 data point 
case. In the other cases the maximum error in the rotation parameter 
estimate is less than nine percent, and for the 100 data point case this 
error is approximately three percent (for the iterative minimization 
scheme) . 
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For higher noise levels (c r = 0. 4 and 0. 5) the one step minimiza- 
tion method runs into serious difficulties. Table 1 shows two instances 
in which the one step minimization method was unable to fit an ellipse 
to the data points. In both of these instances the estimate c? e was found 
to have one of its first two components a negative number. This means 
that the one step minimization method actually fit a hyperbola to the 
given data points rather than an ellipse. 

For relatively high levels of noise (cr = 0. 4 and 0. 5) Table 2 
shows that the iterative minimization scheme also exhibits an undesirable 
characteristic, that is, it has a tendency to select values for the first 
two components of c e which are at the boundary of their respective range 
constraints. When this is the case the resulting fitted ellipse is actually 
a circle having a radius equal to four units, as shown in Figures 12 and 
13. These figures indicate that the iterative minimization scheme has 
attempted to cluster all of the data points along a small portion of the 
boundary of the fitted ellipse ( or circle) , with approximately one half 
of the data points on either side of the fitted ellipse' s boundary. 

4. 5 Summary 

This chapter has investigated the merits of M recognizing n a planar 
elliptical pattern, whose boundary points are given, by estimating the 
values for the five parameters which characterize an ellipse. The 
parameter estimation schemes which were employed are the two which 
were described in Chapter III, namely, the one step minimization method 
and the iterative minimization scheme. The ellipse which was to be 
recognized was permitted to have arbitrary size and shape, as well as 
arbitrary position and orientation so long as it was located within the 
specified field of view. 

As was pointed out in Section 4. 4, the two minimization schemes 
provided essentially the same results for the estimate of the parameter 
vector associated with the reference ellipse when the noise level was 
below cr = 0. 4. For the lower noise levels ( cr = 0. 0, 0.1, and 0.2) these 
estimates were quite good, and special note was made concerning the 
accuracy with which the rotation parameter was estimated. Excluding the 
10 data point case for cr = 0.2, the rotation parameter was never more 
than three percent in error, which is a remarkable result. Unfortunately, 
no other schemes exist presently with which these results can be compared. 

For cr = 0. 3 the estimate for the rotation parameter was still 
respectable, but the other parameters were not estimated accurately 
enough to yield a fitted ellipse which approximated the reference ellipse 
to an acceptable degree. 
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Both minimization schemes displayed undesirable characteristics 
for very high noise levels ( c r =0.4 and 0. 5) . The one step minimization 
method had a tendency to fit a hyperbola to the data points rather than 
an ellipse (characterized by a negative value for one of the first two 
components of the parameter vector) while the iterative minimization 
scheme had a tendency to fit a constrained ellipse to the data points 
( r w = r z = 4. 0) . 

The fact that both minimization schemes failed to accurately 
estimate the parameters associated with the reference ellipse for high 
noise levels does not distract from their usefulness. In practice one 
would regard a system corresponding to cr = 0. 3, 0.4 and 0. 5 as having 
an unacceptable level of noise and hence would demand a better design 
for the system. Upon viewing Figures 5, 6, and 7 one sees that it would 
be very difficult, if not impossible, to develop a recognition scheme that 
could accurately recognize an ellipse from the given scatter of data, 
points (this includes a human being as a "pattern recognizer"). 

As a minor point, it should be mentioned that the undesirable 
characteristics of the two minimization schemes (for high noise) which 
were previously mentioned can be corrected to some extent. The itera- 
tive minimization scheme can be improved if the range constraints on 
the first two components of c^ are further restricted after the data points 
become known. One simple procedure is to construct a rectangle, having 
sides parallel to the x, y-axes, that encloses all the data points and that 
has at least one data point lying on each of its sides. One would expect 
the fitted ellipse to have neither of its diameters larger than the diagonal 
of this "bounding" rectangle. Thus, the two radii are constrained to be 
no larger than one half of this diagonal, and so the lower bounds on the 
range constraints for the first two components of c* e are modified accord- 
ingly. 


For the one step minimization method, constraints could be speci- 
fied so that the fitted pattern is forced to be an ellipse. However, the 
simplicity of the one step minimization method involved the fact that it 
was an unconstrained minimization scheme. Since the unconstrained 
one step minimization method was unable to always fit an ellipse to the 
data points (and for reasons given in Chapter V) the iterative minimiza- 
tion scheme was considered the better scheme and was used for the 
recognition of elliptical patterns as well as for patterns that are not 
ellipses. 


Another point which might be noted concerning the one step 
minimization method is that this scheme tends to estimate the larger 
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radius, r x , much less accurately than the smaller radius, ry, This 
can be seen in Table 1, This is not a property of the iterative minimi- 
zation scheme, however, giving more support for its use. 
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TABLE 1: ELLIPSE PARAMETER ESTIMATES OBTAINED BY 

ONE STEP MINIMIZATION METHOD 


r x = 2. 0 r.y = 1. 0 


O’ 

_N_ 

e i 

® 2 

_A 

_B_ 

e_ 

r x 

r 

_JL 

0. 0 

10 

0. 250035 

1. 000063 

0. 999918 

-1. 999968 

0. 499943 

2. 000 

1. 000 

0. 0 

20 

0. 249987 

0. 999831 

0. 999987 

-2. 000008 

0. 499964 

2. 000 

1. 000 

0. 0 

50 

0. 249962 

1. 000428 

1. 000113 

-2. 000046 

0. 500288 

2. 000 

1. 000 

0. 0 

100 

0. 250081 

1. 000146 

0. 999886 

-1. 999961 

0. 499766 

2. 000 

1. 000 

0. 1 

10 

0. 226588 

0. 985690 

1. 13 3 521 

-2. 017937 

0. 505674 

2. 101 

1.007 

0. 1 

20 

0. 247580 

1. 019186 

0. 986980 

-2. 036948 

0. 500796 

2. 010 

0.991 

0. 1 

50 

0. 230059 

0 . 966594 

1. 061745 

-2. 000900 

0. 545271 

2. 085 

1.017 

0. 1 

100 

0. 240472 

1. 021042 

0. 960186 

-1. 995705 

0. 504009 

2.039 

0.990 

0. 2 

10 

0. 250661 

1. 091854 

1. 046996 

-1. 926707 

0. 610062 

1. 997 

0. 957 

0. 2 

20 

0. 20 5477 

0. 870532 

1. 019036 

-2. 167258 

0. 457618 

2. 206 

1. 072 

0. 2 

50 

0. 245096 

1. 081746 

1. 029216 

-2. 026368 

0. 482127 

2. 020 

0.961 

0. 2 

100 

0. 222750 

1. 104007 

0. 993413 

-2. 054868 

0. 527154 

2. 119 

0. 952 

0. 3 

10 

0. 122579 

1. 141891 

0. 510788 

-2. 504557 

0. 728425 

2. 856 

0. 936 

0. 3 

20 

0. 241930 

1. 038012 

1. 264826 

-1. 812599 

0. 547318 

2. 033 

0. 982 

0. 3 

50 

0. 187579 

0 . 965065 

1. 005195 

-2. 067051 

0. 611433 

2. 309 

1. 018 

0. 3 

100 

0. 141655 

0. 976113 

0. 911619 

-2. 123045 

0. 518779 

2. 657 

1.012 

0. 4 

10 

0. 156311 

0. 662240 

0. 965569 

-1. 958783 

0. 386650 

2. 529 

1. 229 

0. 4 

20 

0. 037682 

1. 751601 

0. 351396 

-2. 735868 

0. 661454 

5. 151 

0. 756 

0. 4 

50 

o. 029652 

0. 997244 

3. 404136 

-1. 132561 

0. 461635 

5. 807 

1. 001 

0. 4 

100 

* 







0. 5 

10 

0. 081115 

0. 644221 

2. 139103 

-1. 733687 

0. 542781 

3. 511 

1. 246 

0. 5 

20 

0. 137709 

0. 967893 

0. 618680 

-2. 566473 

0. 499484 

2. 695 

1.016 

0. 5 

50 

* 







0. 5 

100 

0. 075127 

1. 015188 

1. 068787 

-2. 383365 

0. 507140 

3. 648 

0.992 


’' Unable to fit ellipse. 


Note: Reference ellipse parameters are: e x = 0. 25, e 2 = 1. 00, A = 1. 00, B = -2. 00, 0 = 0. 50 





TABLE 2: 

ELLIPSE PARAMETER 

ESTIMATES OBTAINED BY 






ITERATIVE MINIMIZATION SCHEME 






r x = 2.0 

r y = 1.0 






(T 

N 

e t 

e 2 

A 

B 

e 

r T 

r 

* 










0. 0 

10 

. 25000 

1.00000 

1.00000 

-2. 00000 

0. 50000 

2. 000 

1. 000 

0. 0 

20 

. 25000 

1. 00000 

1. 00000 

-2. 00000 

0. 50000 

2. 000 

1. 000 

0. 0 

50 

. 25000 

1. 00000 

1. 00000 

-2. 00000 

0. 50000 

2. 000 

1. 000 

0. 0 

100 

. 25000 

1.00000 

1.00000 

-2. 00000 

0. 50000 

2. 000 

1.000 

0. 1 

10 

. 23309 

0. 91828 

1. 12770 

-2. 00796 

0. 50372 

2. 071 

1. 044 

0. 1 

20 

. 25481 

0. 97564 

0. 98513 

-2. 03016 

0. 50099 

1. 981 

1. 012 

0. 1 

50 

. 24045 

0. 90386 

1. 05204 

-1. 99382 

0. 54868 

2. 039 

1. 052 

0. 1 

100 

. 25103 

0. 95177 

0. 95558 

-1. 98589 

0. 50 581 

1.996 

1. 025 

0. 2 

10 

. 26219 

0. 98635 

1. 04461 

-1. 91310 

0. 60783 

1. 953 

1. 007 

0. 2 

20 

. 23523 

0. 65429 

0. 99513 

-2. 13235 

0. 46806 

2. 062 

1. 236 

0. 2 

50 

. 28036 

0. 86192 

1. 02352 

-1. 98737 

0 . 49220 

1. 889 

1. 077 

0. 2 

100 

. 26151 

0. 86500 

0. 97935 

-2. 01866 

0. 54554 

1. 955 

1. 075 

0. 3 

10 

. 17273 

0. 78372 

0. 56804 

-2. 39295 

0. 80087 

2. 406 

1. 130 

0. 3 

20 

. 32107 

0. 62374 

1. 13235 

-1. 81673 

0. 61739 

1. 765 

1. 266 

0. 3 

50 

. 23160 

0. 65681 

0. 99362 

-1. 99064 

0. 63811 

2. 078 

1. 234 

0. 3 

100 

. 22484 

0. 60177 

0. 88051 

-2. 01215 

0. 54895 

2. 109 

1. 289 

0. 4 

10 

. 16435 

0. 53921 

0. 80733 

-1. 93123 

0. 30713 

2. 467 

1. 362 

0. 4 

20 

. 20117 

0. 74273 

0. 79921 

-2. 26838 

0. 70981 

2. 230 

1. 160 

0. 4 

50 

. 06250 

0. 06250 

-0. 55830 

1. 00064 

0. 11353 

4. 000 

4. 000 

0. 4 

100 

. 06250 

0. 06250 

-1. 12144 

1. 08820 

0. 12380 

4. 000 

4. 000 

0. 5 

10 

. 10609 

0. 55384 

1. 85071 

-1. 87407 

0. 49956 

3.070 

1. 344 

0. 5 

20 

. 26380 

0. 44969 

0. 74672 

-2. 38593 

0. 64307 

1. 947 

1.491 

0. 5 

50 

. 06250 

0. 06250 

-0. 04203 

1. 33321 

0. 13402 

4. 000 

4. 000 

0. 5 

100 

. 06250 

0. 06250 

-1. 35350 

0. 90174 

0. 14893 

4. 000 

4. 000 


*Note: 

Reference ellipse parameters are: e x = 

0. 25, e 2 = 1. 00, 

A = 1. 00, B = 

2.00, e 

= 0. 50 
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CHAPTER V 


THE RECOGNITION OF RECTANGULAR PLANAR PATTERNS 
5. 1 Introduction 

The objective of this chapter is to investigate the feasibility of 
employing either or both of the estimation schemes which are discussed 
in Chapter III to recognize rectangular planar patterns. Again, by the 
term "recognition 11 is meant the estimation of the five parameters which 
characterize a rectangle having arbitrary size and shape, as well as 
arbitrary position (translation and rotation) in a planar field of view. 

In Section 5. 2 the statement' of the problem is formulated and the 
recognition strategy which will be investigated is discussed. The results 
which are obtained from the two minimization schemes for rectangles 
having known parameter vectors (no noise) are then analyzed in Section 
5. 3. 


The iterative minimization scheme is employed in Section 5.4 to 
estimate the parameters associated with rectangles whose boundary 
points are corrupted with noise, and the results using this scheme are 
discussed. Finally, Section 5. 5 contains a brief summary of the results 
and conclusions which have been reached in this chapter. 

5. 2 Statement of Problem 

In order to further test the recognition techniques developed for 
elliptical objects, a second class of patterns was considered. Rectangu- 
lar patterns were chosen for this purpose because they are simple geo- 
metric patterns and yet do not possess a simple analytic representation. 
Furthermore, a rectangle has several properties in common with an 
ellipse. Both of these patterns are convex, and both are symmetrical 
about two orthogonal axes. Thus, a rectangle may be characterized by 
a set of five parameters in a manner quite similar to an ellipse. Figure 
14 shows a rectangle which has been translated and rotated with respect 
to the x-y coordinate system. This rectangle is characterized by its 
w-axis radius R w and z-axis radius R z , and by the x- and y-translation of 
its center (A 1 and B f , respectively) , as well as by its rotation 0 1 . 
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Therefore, the parameter vector characterizing a rectangle may be 
expressed as 



B' 

6 ' 


( 5-1) 


The strategy which was used to recognize rectangular patterns 
was to present a number of different size and shape rectangles to both 
the one step minimization method and the iterative minimization scheme, 
and to determine what relationship, if any, existed between the parameter 
vectors of the fitted ellipses and the parameter vectors of the correspond- 
ing rectangles. This strategy is motivated by the fact that the noiseless 
data points which lie on the boundary of a rectangle may be considered as 
being noisy data points which originally belonged on the boundary of some 
ellipse. If a relationship can be found between the parameter vectors of 
the fitted ellipses and the parameter vectors of the corresponding rec- 
tangles, then it will be possible to compute the parameter vector of an 
unknown rectangle after the parameter vector of its associated fitted 
ellipse is determined. 


5. 3 Parameter Estimation for Noise- 
Free Rectangular Patterns 

The results of using the one step minimization method to recognize 
a rectangle are shown in Tables 3, 4, 5, and 6. Each table corresponds 
to a different number of data points on the boundary of the rectangle, 
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the number of data points being 8, 20, 48, or 100, respectively. These 
data points were generated by the DATA subroutine which is described 
in Appendix I. 

A total of twenty rectangles were to be recognized. Ten of the 
rectangles have a w-axis radius equal to one unit of length, with the 
z-axis radius varying from 0. 1 to 1.0 units of length in increments of 0. 1. 
The other ten rectangles are exactly twice the dimensions of the first 
ten. The translation and rotation parameters of the rectangles were 
chosen to be the same as those which were used for the ellipse which 
was discussed in Chapter IV. Thus 

Cj = x-translation = A 1 =1.0 
c i = y-translation = B 1 = - 2. 0 

c^ = rotation in radians = 0' = 0. 5 ( 5. 2) 

In the tables the radii of the rectangles are denoted by R while 
the radii of the fitted ellipses are denoted by r. The radii of the fitted 
ellipses are, of course, related to the first two components of their 
characterizing parameter vector, c?, by Eqn. ( 5. 3) . 



( 5.3) 

Inspection of Tables 3, 4, 5, and 6 indicates that the one step 
minimization method is not very effective in recognizing the rectangles. 
For the most part the translation and rotation parameters of the fitted 
ellipses are quite close in value to the corresponding parameters of the 
given rectangles. However, there is no recognizable correspondence 
between the w-axis and z-axis radii of the fitted ellipses and the respec- 
tive radii of the rectangles. 

A remark should be made at this point regarding the results 
which one would intuitively expect to obtain. First of all, since an ellipse 
and rectangle have similar geometric properties ( convexity and symmetry 
abom two orthogonal axes) one would expect that the known rectangles 
and the corresponding fitted ellipses would have identical coordinates for 
their centers, as well as identical rotation angles. On the other hand, 
it is difficult to predict the exact relation between each radius of a known 
rectangle and the corresponding radius of the fitted ellipse. However, 
again due to symmetry, one would expect that the ratio of the radii of a 
fitted ellipse would be nearly equal to the ratio of the radii of the cor- 
responding known rectangle, being more or less independent of the rec- 
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tangle' s size. The one step minimization method did not possess this 
property, however. If a recognition scheme does have this property, 
then the constant of proportionality relating the size of the rectangle to 
the size of the corresponding fitted ellipse may be determined experi- 
mentally. 


Another weakness of the one step minimization method is that it 
is unable to fit any ellipse to some of the given rectangles. This might 
be expected, however, since the one step minimization was unable to 
fit an ellipse to an ellipse under high noise conditions, as was pointed 
out in Chapter IV. Thus, the one step minimization method does not 
appear to be a good method for estimating the parameters of a rectangle. 


The iterative minimization scheme was next employed to esti- 
mate the parameters of these same rectangles. The same values were 
used for the range constraints as were used in the recognition of ellipses 
in Chapter IV. Tables 7, 8, 9, and 10 show the results of using this 
scheme. The initial estimate used for the parameter vector of the fitted 
ellipse for the ten larger rectangles was 


"0. 25*" 
1. 00 
0. 7 
- 2. 5 
. 0.8 _ 


(5.4) 


Inspection of Tables 7, 8, 9 , and 10 indicates that the iterative 
minimization scheme is quite effective in estimating the parameters of 
a rectangle. It is seen that the ratio r z /r w corresponding to the radii 
of the fitted ellipse is equal (to within three decimal places) to the radii 
R 2 /R w rectangle which is to be recognized. Thus the ratio of 

the radii of the fitted ellipses gives a direct indication of the shape of 
the rectangles to which they are fitted. 

Tables 7, 8, 9> and 10 also indicate that the size of the rec- 
tangles may be determined to a reasonable degree of accuracy. When 
eight data points on the boundary of the rectangle are given, Table 7 
shows that R w = 0. 77 5 r w , meaning that the rectangles' radii are 0. 775 
times the length of the fitted ellipses' radii. Likewise, Tables 8, 9 , 
and 10 show that the scale factor, R w /r w is equal to 0. 830, 0. 842, and 
0. 845 for 20, 48, and 100 data points on the boundary of the rectangle, 
respectively. 

One can see that the scale factor does not change appreciably 
when more than 48 data points are given. If one assumes that the scale 
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factor associated with 100 data points is essentially the same as the 
scale factor associated with an infinite number of data points (which seems 
reasonable in light of the above results) then it is possible to compare 
the scale factor associated with a finite number of data points with the 
scale factor associated with a continuous representation of the rectangle. 
For eight data points this ratio is 0. 775/0. 845 = 0. 917, which means 
that the estimated size of the rectangle is only 91- 7% of the size of the 
actual rectangle, although it has exactly the same shape as the actual 
rectangle. For 20 data points this ratio increases to 0. 830/0. 845 = 

0. 983 and for 48 data points the ratio is 0. 842/0. 845 = 0, 997. Thus, if 
the rectangle is represented by twenty or more data points on its boundary, 
then one need merely multiply the radii of the fitted ellipse by the factor 
0. 845 to obtain the radii of the corresponding rectangle, having assurance 
that this rectangle will be at least within 2% of the size of the actual 
rectangle. 

Tables 7, 8, 9? and 10 show another very desirable property of 
the iterative minimization scheme, namely, the translation and rotation 
parameters of the fitted ellipses have exactly the same values (to within 
three decimal places) as the corresponding parameters of the given 
rectangles. Therefore, only the first two components of the fitted 
ellipse 1 s parameter vector need to be transformed in order to obtain the 
parameter vector of the rectangle, and this transformation is a simple 
scale change. 

Thus, the desired relationship between the parameter vector of 
the fitted ellipse and the parameter vector of the associated rectangle 
has now been determined. If the final estimate for the parameter vector 
of the fitted ellipse is given by 



( 5. 5) 


then the estimate for the parameter vector of the rectangle which is to 
be recognized is 


R w 


k/\T 

R z 


k / nT e 2 

A’ 

= 

A 

B' 


B 

0' . 


.0 
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where k, the scale factor, is a function of the number of data points. 

If the number of data points is 20 or greater, k may be taken to be 0. 845. 

5. 4 Parameter Estimation for Noisy 
Rectangular Patterns 

Since the iterative minimization scheme was able to effectively recognize 
rectangles, the question naturally arises as to how well it can recognize 
rectangles which are represented by noisy data points. In order to deter- 
mine this a total of twenty-four different cases were considered, as was 
done with the ellipse in Chapter IV. 

The noiseless rectangle, whose parameters are to be estimated, 
is characterized by the following parameter vector. 


B w 


r 2. 00‘ 

*z 


1. 00 

A' 

= 

1. 00 

B 1 


-2. 00 

_ 0 ' . 


_0. 50. 


The last three components have the same value as they did for the 
ellipse considered in Chapter IV. The noise levels are also the same as 
they were previously, namely, o* = 0. 0, 0. 1, 0. 2, 0. 3, 0. 4, and 0. 5. Also, 
all of the parameters associated with the iterative minimization scheme 
were given the same values as they had in Chapter IV, with the exception 
of the initial parameter vector estimate <? e . It is 


e i 


■ 0. 25] 

e 2 


1. 00 

A 

= 

O 

O 

B 


-2. 50 

.0 - 


o 

00 

o 


Figures 15, 16, 17, 18, 19, and 20 show the results of estimating 
the parameter vector of a rectangle using the iterative minimization 
scheme. Before discussing these results it should be pointed out that 
the radii of the rectangle were computed by using the scale factor associ- 
ated with the appropriate number of data points. Thus, for example, for 
the six cases in which the rectangle was represented by 20 data points 
the scale factor which was used was 0. 830, and not 0. 845. By doing this, 
any error in the estimated parameter vector is due to the noisy data points. 


An examination of Figure 2 reveals that the x-radii and y-radii 
(w-axis radii and x-axis radii, respectively) of the fitted rectangles 
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differ in the third decimal place from the corresponding radii of the 
reference rectangles. This is due to rounding off the scale factors to 
the third decimal place. This error is entirely negligible compared to 
the error resulting from the noisy data points. 

In some of the figures there are not as many noisy data points 
as the number which is indicated. This is due to the fact that some of 
the noisy data points fell outside of the field of view. As before, these 
data points are considered to be valid points for the iterative minimiza- 
tion scheme to use, being analogous to measurement noise. 

An inspection of Figures 15, 16, 17, 18, 19, and 20 shows that 
fitting an ellipse to a set of data points belonging on the boundary of a 
rectangle is an effective method by which to estimate the parameters 
of the rectangle when the noise level is within reasonable limits 
(cr =0.0, 0.1, and 0.2). 

Referring to Figure 16 ( or = 0. 1) one sees that the rectangles 
which correspond to the fitted ellipses (the dashed line rectangles) 
resemble the reference rectangles very closely except in the eight data 
point case. It would seem that eight data points, when corrupted by 
noise, simply are too sparse in number for the iterative minimization 
scheme to yield a good estimate for the reference rectangle 1 s parameter 
vector. However, it should be noted that the error in estimating the 
rotation angle for the eight data point case is quite acceptable, being 
approximately two percent. In the other three cases this error is approxi- 
mately one percent or less. 

For cr = 0. 2 it can be seen in Figure 17 that the parameter vector 
estimates are beginning to deteriorate, but for the 20, 48 and 100 data 
point cases these estimates are still acceptable by most standards. In 
particular, it is seen that for these three cases the error in the rotation 
angle estimate is no greater than approximately four percent, which is 
rather small considering the scatter of the data points. 

When the noise level reaches cr =0.3, the overall effectiveness 
of the iterative minimization scheme becomes questionable. The fitted 
rectangles have a tendency to be larger than the reference rectangles. 
However, one good point which can be made is that the estimate for the 
rotation angle in all four cases does not exceed four percent, which 
means that this estimate has not been affected to any extent by the 
increase in noise level from ( r = 0. Z to cr = 0. 3. 
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Figures 19 and 20 indicate that for high noise levels ( cr s 0.4 and 
0. 5) the recognition capability of the iterative minimization scheme has 
completely deteriorated. Some improvement could be achieved for the 
cases in which the fitted rectangles have radii equal to their constraint 
value. In these cases the fitted rectangle is a square wm.cn concentrates 
the data points in one of its corners, witn approximately one half of the 
data points on ttie inside of the square and one half on the outside. This 
situation is very similar to that which occurred in the recognition of 
ellipses under high noise conditions, and it can be remedied in exactly 
the same manner as described in Chapter IV. 

5. 5 Summary 

This chapter investigated the feasibility of utilizing either the 
one step minimization method or the iterative minimization scheme to 
estimate the parameters of a rectangle when noisefree data points lying 
on the rectangle 1 s boundary are given. If the parameters are estimated 
with sufficient precision then the rectangle has been M recognized M correctly. 

It was found that the one step minimization method was completely 
inadequate in its capability to estimate the parameters of given rectangles. 
While it did do a reasonable job in estimating the translation and rotation 
parameters, the two major shortcomings of this method were 

( 1) the ellipse which was fitted to the data points did not 
have the same shape as the given rectangle, i. e. , the 
ratio of the radii of the fitted ellipse was not identical 
to the corresponding ratio of the radii of the given 
rectangle , 

and 

( 2) the size of the fitted ellipse did not double when the size of 
the given rectangle doubled. 

The iterative minimization scheme, on the other hand, did not 
have these shortcomings. Not only did it estimate the translation and 
rotation parameters very precisely, but the ellipse which it fitted to the 
data points had radii whose ratio was identical to that of the given 
rectangle, and this ratio was .independent of the size of the given rec- 
tangle. It was therefore possible to experimentally determine a scale 
factor relating the size of the fitted ellipse 1 s radii to the radii of the 
given rectangle. 

Since the iterative minimization scheme had the capability to 
precisely estimate the parameters of a rectangle whose boundary points 
were noise free, the next step was to determine the degradation in the 
parameter vector estimates in situations for which the data points were 
noisy. Reference to Figures 15, l6, and 17 indicates that for moderate 
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levels of noise ( <r = 0. 0, 0. 1, and 0. 2) the iterative minimization scheme 
did a very satisfactory job of recognizing the rectangles. Special notice 
should be taken concerning the accuracy with which the rotation angle 
was estimated. Excluding the eight data point case, this error was 
never greater than four percent for these moderate noise levels. 

The recognition scheme produced results of questionable value 
for noise level o' = 0. 3. Although the rotation angle was still estimated 
with good precision (maximum of four percent error) , the size of the 
fitted rectangle tended to be larger than the size of the reference rec- 
tangle. It can be said that <r = 0. 3 represents the maximum noise level 
for which the iterative minimization scheme produces useful results for 
the particular set of rectangles investigated. 

For larger noise levels (cr = 0.4 and 0. 5) the iterative minimiza- 
tion scheme was not able to do a satisfactory job of recognizing the rec- 
tangles at all. This is not at all surprising, since even a human being 
would have difficulty trying to fit a rectangle to the data as shown on 
Figures 19 and 20. 
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TABLE 3: PARAMETERS OF ELLIPSES FITTED TO RECTANGLES 

BY THE ONE STEP MINIMIZATION METHOD 
(8 DATA POINTS) 




® Z^W 

r w 

r z 

r z^ r w 

A 

B 

0 

1 . 0 

0 . 1 

0 . 100 

1. 331 

0. 077 

. 058 

1. 152 

-1. 917 

. 496 

1 . 0 

0 . 2 

0 . 200 

1. 768 

0 . 202 

. 114 

. 962 

- 2 . 023 

. 499 

1 . 0 

0. 3 

0. 300 

1. 755 

0 . 298 

. 170 

. 996 

-2. 007 

. 499 

1 . 0 

0. 4 

0. 400 

1. 773 

0 . 396 

. 223 

1. 003 

- 2 . 008 

. 499 

1 . o 

0. 5 

0. 500 

1. 799 

0 . 491 

. 273 

1 . 006 

-2. 016 

. 498 

1 . o 

0 . 6 

0 . 600 

1. 834 

0. 585 

. 319 

1. 007 

- 2 . 018 

. 497 

1 . o 

0. 7 

0. 700 

1 . 886 

0. 675 

. 358 

1 . 010 

-2. 025 

. 496 

1.0 

0 . 8 

0 . 800 

1. 957 

0 . 761 

. 389 

1. 014 

-2. 034 

. 49 5 

1 . 0 

0 . 9 

0 . 900 

2. 063 

0. 842 

. 408 

1.017 

-2. 045 

. 493 

1 . o 

1 . 0 

1 . 000 

2. 234 

0. 915 

. 410 

1 . 020 

- 2 . 060 

. 491 

2 . 0 

0 . 2 

0 . 100 

3. 517 

0 . 201 

. 057 

. 995 

-2. 005 

. 500 

2 . 0 

0. 4 

0 . 200 

3. 540 

0. 395 

. 112 

1. 003 

-2. 007 

. 500 

2 . 0 

0 . 6 

0. 300 

3. 670 

0. 584 

. 160 

1 . 008 

-2. 017 

. 499 

2 . 0 

0 . 8 

0. 400 

3. 918 

0 . 761 

. 194 

1. 014 

-2. 034 

. 499 

2 . 0 

1 . 0 

0. 500 

4. 481 

0. 915 

. 204 

1 . 021 

- 2 . 060 

. 498 

2 . 0 

1 . 2 

0 . 600 

6 . 801 

1. 030 

. 151 

1 . 010 

- 2 . 109 

. 497 

2 . 0 

1. 4 

0. 700 

* 






2 . 0 

1 . 6 

0 . 800 

* 






2 . 0 

1 . 8 

0 . 900 







2 . 0 

2 . 0 

1 . 000 

n 5 







Not able to fit ellipse data. 

Note: R corresponds to reference rectangles' radii 

r corresponds to fitted ellipses' radii 
All reference rectangles have A = 1.0, B = -2.0, and 

0 = 0.5 
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TABLE 4: PARAMETERS OF ELLIPSES FITTED TO RECTANGLES 

BY THE ONE STEP MINIMIZATION METHOD 
(20 DATA POINTS) 



® z 

R z/^ W 

r w 

r z 

r z / r w 

A 

B 

0 

1. 0 

0. 1 

0. 100 

0. 699 

0. 056 

0. 080 

1. 238 

-1. 865 

0. 498 

1. o 

0. 2 

0. 200 

1. 326 

0. 209 

0. 158 

0. 988 

-2. 007 

0. 500 

1. o 

0. 3 

0. 300 

1. 330 

0. 314 

0. 236 

1. 003 

-2. 002 

0. 500 

1. 0 

0. 4 

0. 400 

1. 335 

0. 418 

0. 313 

1. 003 

-2. 004 

0. 499 

1. o 

0. 5 

0. 500 

1. 340 

0. 520 

0. 388 

1. 004 

-2. 008 

0. 498 

1. o 

0. 6 

0. 600 

1. 347 

0. 622 

0. 462 

1. 005 

-2. 012 

0. 497 

1. 0 

0. 7 

0. 700 

1. 354 

0. 722 

0. 533 

1. 008 

-2. 016 

0 . 496 

1. 0 

0. 8 

0. 800 

1. 364 

0. 819 

0. 600 

1. Oil 

-2. 022 

0. 494 

1.0 

0. 9 

0. 900 

1. 377 

0. 914 

0. 664 

1. 014 

-2.029 

0. 491 

1. o 

1.0 

1. 000 

1. 394 

1. 005 

0. 721 

1. 018 

-2. 038 

0. 487 

2. 0 

0. 2 

0. 100 

2. 643 

0. 209 

0. 079 

1. 017 

-1. 992 

0. 500 

2. 0 

0. 4 

0. 200 

2. 671 

0. 418 

0. 156 

1. 002 

-2. 005 

0. 500 

2. 0 

0. 6 

0. 300 

2. 69 5 

0. 622 

0. 231 

1. 005 

-2. 012 

0. 499 

2. 0 

0. 8 

0. 400 

2. 729 

0. 819 

0. 300 

1. 011 

-2. 022 

0. 499 

2. 0 

1. o 

0. 500 

2. 789 

1. 004 

0. 360 

1. 018 

-2. 038 

0. 498 

2. 0 

1. 2 

0. 600 

2. 902 

1. 168 

0. 402 

1. 029 

-2. 063 

0. 497 

2. 0 

1. 4 

0. 700 

3. 165 

1. 289 

0. 407 

1. 046 

-2. 104 

0. 49 5 

2. 0 

1. 6 

0. 800 

4. 429 

1. 309 

0. 296 

1. 062 

-2. 187 

0 . 492 

2. 0 

1. 8 

0. 900 

* 






2. 0 

2. 0 

1. 000 

£ 







n Not able to fit ellipse to data. 
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TABLE 5: PARAMETERS OF ELLIPSES FITTED TO RECTANGLES 

BY THE ONE STEP MINIMIZATION METHOD 
(4 8 DATA POINTS) 


R w 

Rz 

® z 

r w 

r z 

r z/r w 

A 

B 

e 

1. 0 

0. 1 

0. 100 

1. 307 

0. Ill 

0. 085 

0. 640 

-2. 206 

0. 515 

1. 0 

0. 2 

0. 200 

1. 246 

0. 207 

0. 166 

1. 001 

-2. 001 

0. 503 

1. 0 

0. 3 

0. 300 

1. 277 

0. 318 

0. 249 

1. 001 

-2. 002 

0. 501 

1. 0 

0. 4 

0. 400 

1. 282 

0. 424 

0. 331 

1. 001 

-2. 005 

0. 500 

1.0 

0. 5 

0. 500 

1. 284 

0. 529 

0. 412 

1. 003 

-2. 007 

0 . 499 

1.0 

0. 6 

0. 600 

1. 288 

0. 632 

0 . 491 

1. 004 

-2. 010 

0. 498 

1. o 

0. 7 

0. 700 

1. 292 

0. 735 

0 . 569 

1. 006 

-2. 013 

0. 497 

1. 0 

0. 8 

0. 800 

1. 299 

0. 835 

0. 643 

1. 009 

-2. 018 

0. 495 

1.0 

0. 9 

0. 900 

1. 307 

0. 933 

0. 714 

1. 012 

-2. 024 

0 . 491 

1.0 

1. o 

1. 000 

1. 317 

1. 028 

0. 781 

1.015 

-2. 031 

0. 486 

2. 0 

0. 2 

0. 100 

2. 537 

0. 21 1 

0. 083 

0. 977 

-2. 014 

0. 500 

2. 0 

0. 4 

0. 200 

2. 564 

0. 424 

0. 165 

1 . 002 

-2. 004 

0. 500 

2. 0 

0. 6 

0. 300 

2.576 

0. 633 

0. 246 

1. 004 

-2. 010 

0. 500 

2. 0 

0. 8 

0. 400 

2. 598 

0. 835 

0. 321 

1. 009 

-2. 018 

0. 499 

2. 0 

1. o 

0. 500 

2. 635 

1 . 028 

0. 390 

1. 015 

-2. 032 

0. 498 

2. 0 

1 . 2 

0. 600 

2. 702 

1. 203 

0. 445 

1. 025 

-2. 052 

0. 497 

2. 0 

1. 4 

0. 700 

2. 845 

1. 343 

0. 472 

1. 040 

-2. 086 

0. 495 

2. 0 

1 . 6 

0. 800 

3. 323 

1. 401 

0. 422 

1. 065 

-2. 152 

0 . 492 

2. 0 

1 . 8 

0. 900 

* 






2. 0 

2. 0 

1 . 000 

* 







& 

Not able to fit ellipse to data. 
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TABLE 6 ; PARAMETERS OF ELLIPSES FITTED TO RECTANGLES 
BY THE ONE STEP MINIMIZATION METHOD 
(100 DATA POINTS) 


R w 

Rz 

® z^^w 

r w 

r z 

r z /r w 

A 

B 

e 

1. o 

0. 1 

0. 100 

0. 491 

0. 044 

0. 090 

0. 722 

-2. 141 

0. 506 

1.0 

0. 2 

0. 200 

1. 256 

0. 213 

0. 170 

1.015 

-1. 992 

0. 503 

1. o 

0. 3 

0. 300 

1. 265 

0. 320 

0. 253 

0. 995 

-2. 004 

0. 501 

1. o 

0. 4 

0. 400 

1. 269 

0. 425 

0. 335 

1. 000 

-2. 004 

0. 500 

1.0 

0. 5 

0. 500 

1. 273 

0. 530 

0. 416 

1. 002 

-2. 006 

0. 499 

1. o 

0. 6 

0. 600 

1. 277 

0. 635 

0. 497 

1. 004 

-2. 009 

0. 498 

1. o 

0. 7 

0. 700 

1. 281 

0. 737 

0. 575 

1. 006 

-2. 013 

0 . 496 

1. 0 

0. 8 

0. 800 

1. 28 7 

0. 838 

0. 651 

1. 008 

-2. 017 

0. 495 

1. o 

0. 9 

0. 900 

1. 29 5 

0. 937 

0. 724 

1.011 

-2. 023 

0 . 491 

1. 0 

1. 0 

1. 000 

1. 305 

1. 032 

0 . 791 

1. 015 

-2. 030 

0. 486 


2. 0 

0. 2 

0. 100 

2. 552 

0. 214 

0. 084 

0. 995 

-2. 005 

0. 497 

2. 0 

0. 4 

0. 200 

2. 543 

0. 425 

0. 167 

1. 003 

-2. 004 

0. 499 

2. 0 

0. 6 

0. 300 

2. 556 

0. 634 

0. 248 

1. 005 

-2. 009 

0 . 499 

2. 0 

0. 8 

0. 400 

2. 577 

0. 838 

0. 325 

1. 008 

-2. 018 

0 . 499 

2. 0 

1. 0 

0. 500 

2. 611 

1. 032 

0. 395 

1. 015 

-2. 030 

0. 498 

2. 0 

1. 2 

0. 600 

2. 672 

1. 209 

0. 452 

1. 024 

-2. 050 

0. 497 

2. 0 

1. 4 

0. 700 

2. 801 

1. 353 

0. 483 

1. 039 

-2. 083 

0. 495 

2. 0 

1. 6 

0. 800 

3. 215 

1. 417 

0. 441 

1. 064 

-2. 146 

0 . 492 

2.0 

1. 8 

0 . 900 

-r* 






2. 0 

2. 0 

1. 000 








n Not able to fit ellipse to data. 
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TABLE 7: PARAMETERS OF ELLIPSES FITTED TO RECTANGLES 

BY THE ITERATIVE MINIMIZATION METHOD 
(8 DATA POINTS) 



R z 

F z / R w 

r w 

r z 

r z/ r w 

A 

B 

e 

2. 0 

0. 2 

0. 100 

2. 582 

0. 258 

0. 100 

1. 000 

- 2. 000 

0. 500 

2.0 

0. 4 

0. 200 

2. 582 

0. 516 

0. 200 

1. 000 

- 2. 000 

0. 500 

2. 0 

0. 6 

0. 300 

2. 582 

0. 775 

0. 300 

1. 000 

- 2, 000 

0. 500 

2. 0 

0. 8 

0. 400 

2. 582 

1. 033 

0. 400 

1. 000 

-2. 000 

0. 500 

2. 0 

1. o 

0. 500 

2 . 582 

1. 291 

0. 500 

1. 000 

-2. 000 

0. 500 

2. 0 

1. 2 

0. 600 

2. 582 

1. 549 

0. 600 

1. 000 

-2. 000 

0. 500 

2. 0 

1. 4 

0. 700 

2. 58 2 

1. 807 

0. 700 

1. 000 

-2. 000 

0. 500 

2. 0 

1. 6 

0. 800 

2. 582 

2. 066 

0. 800 

1. 000 

-2. 000 

0. 500 

2. 0 

1. 8 

0. 900 

2. 582 

2. 324 

0. 900 

1. 000 

- 2. 000 

0. 500 

2. 0 

2. 0 

1. 000 

2. 582 

2. 582 

1. 000 

1. 000 

- 2. 000 

0. 454 


r - ellipse radii 
R - rectangle radii 


Scale Factor 


= k 8 


|8 data points 



r w r z 


0. 775 
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TABLE 8: PAEAMETEES OF ELLIPSES FITTED TO RECTANGLES 

BY THE ITERATIVE MINIMIZATION METHOD 
( 20 DATA POINTS) 


R w 


R 

r w 

r z 

r z /r w 

A 

B 

0 

2. 0 

0. 2 

0. 100 

2. 409 

0. 241 

0. 100 

1. 000 

-2. 000 

0. 500 

2. 0 

0. 4 

0. 200 

2. 409 

0. 482 

0. 200 

1. 000 

-2. 000 

0. 500 

2. 0 

0. 6 

0. 300 

2. 409 

0. 723 

0. 300 

1. 000 

-2. 000 

0. 500 

2. 0 

0. 8 

0. 400 

2. 409 

0. 963 

0. 400 

1. 000 

-2. 000 

0. 500 

2. 0 

I. o 

0. 500 

2. 409 

1. 204 

0. 500 

1. 000 

-2. 000 

0. 500 

2. 0 

1. 2 

0. 600 

2. 409 

1. 445 

0. 600 

1. 000 

-2. 000 

0. 500 

2. 0 

1. 4 

0. 700 

2. 409 

1. 686 

0. 700 

1. 000 

-2. 000 

0. 500 

2, 0 

1. 6 

0. 800 

2. 409 

1. 927 

0. 800 

1. 000 

-2. 000 

0. 500 

2. 0 

1. 8 

0. 900 

2. 409 

2. 168 

0. 900 

1. 000 

-2. 000 

0. 500 

2. 0 

2. 0 

1. 000 

2. 409 

2. 409 

1. 000 

1. 000 

-2. 000 

0. 859 


r - ellipse radii 


R - rectangle radii 

= k 20 = Rw = Pz =0. 830 

Scale F actor! 


r w r z 

20 data points 
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TABLE 9: PARAMETERS OF ELLIPSES FITTED TO RECTANGLES 

BY THE ITERATIVE MINIMIZATION METHOD 
( 48 DATA POINTS) 


R w 

Rz 

^ z w 

r w 

r z 

r z/ r w 

A 

B 

0 

2. 0 

0. 2 

0. 100 

2. 374 

0. 237 

0. 100 

1. 000 

-2. 000 

0, 500 

2. 0 

0. 4 

0, 200 

2. 374 

0. 475 

0. 200 

1. 000 

- 2, 000 

0. 500 

2. 0 

0, 6 

0, 300 

2. 374 

0. 712 

0. 300 

1. 000 

-2. 000 

0. 500 

2. 0 

0. 8 

0. 400 

2, 374 

0. 950 

0, 400 

1. 000 

-2. 000 

0, 500 

2. 0 

1. 0 

0. 500 

2. 374 

1. 187 

0. 500 

1, 000 

-2. 000 

0. 500 

2. 0 

1. 2 

0. 600 

2. 374 

1. 424 

0. 600 

1. 000 

-2. 000 

0. 500 

2. 0 

1. 4 

0. 700 

2. 374 

1. 662 

0. 700 

1. 000 

-2, 000 

0. 500 

2. 0 

1. 6 

0. 800 

2. 374 

1. 899 

0. 800 

1. 000 

-2. 000 

0. 500 

2. 0 

1. 8 

0. 900 

2. 374 

2. 137 

0. 900 

1. 000 

-2, 000 

0. 500 

2. 0 

2. 0 

1. 000 

2. 374 

2. 374 

1. 000 

1. 000 

-2. 000 

0. 864 


r - ellipse radii 
R - rectangle radii 
Scale Factorl 


= k 


R w R z 


48 " 


= 0, 842 


r r 
w z 


48 data points 
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TABLE 10; PARAMETERS OF ELLIPSES FITTED TO RECTANGLES 
BY THE ITERATIVE MINIMIZATION METHOD 
(100 DATA POINTS) 


R w 

R z 


r w 

r z 

r z/ r W 

A 

B 

0 

2. 0 

0. 2 

0. 100 

2. 368 

0. 237 

0. 100 

1. 000 

-2. 000 

0. 500 

2. 0 

0. 4 

0. 200 

2. 368 

0. 474 

0. 200 

1. 000 

-2. 000 

0, 500 

2. 0 

0. 6 

0. 300 

2. 368 

0. 710 

0. 300 

1. 000 

-2. 000 

0. 500 

2, 0 

0 . 8 

0. 400 

2. 368 

0. 947 

0. 400 

1. 000 

-2, 000 

0. 500 

2. 0 

1. o 

0. 500 

2. 368 

1. 184 

0, 500 

1. 000 

-2. 000 

0. 500 

2. 0 

1. 2 

0. 600 

2. 368 

1. 421 

0. 600 

1. 000 

-2. 000 

0. 500 

2. 0 

1. 4 

0. 700 

2. 368 

1. 658 

0. 700 

1. 000 

-2. 000 

0. 500 

2. 0 

1. 6 

0. 800 

2. 369 

1. 895 

0. 800 

1. 000 

-2. 000 

0. 500 

2. 0 

1. 8 

0. 900 

2. 368 

2. 131 

0. 900 

1. 000 

-2. 000 

0 . 500 

2, 0 

2. 0 

1. 000 

2, 368 

2. 368 

1. 000 

1. 000 

-2. 000 

0, 500 


r - ellipse radii 


R - rectangle radii 
Scale Factor = k x 0 o 

100 data points 



0. 845 
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CHAPTER VI 


SUMMARY AND CONCLUSIONS 


This investigation has been aimed at a solution of the problem of 
real time landmark identification from spacecraft optical fields. The 
approach which has been taken relies upon the reduction of two dimen- 
sional optical images to a discrete set of data points associated with the 
boundary of an object to be identified. Granted such a set of points, the 
work reported here is directed toward the fitting of computationally 
generated images to the real image points by means of an algorithm based 
upon nonlinear regression analysis. 

In order to obtain some concrete results, the present investigation 
has been limited to a consideration of arbitrary elliptical and arbitrary 
rectangular objects. While such objects may be relatively rare among 
all possible landmarks of interest, the approach taken is one of approxi- 
mation of irregular objects by elliptical or rectangular templates. That 
is, the methods developed are tolerant of large amounts of noise whether 
this noise is introduced by measurement and sensing processes or by 
the deviation of real objects from elliptical or rectangular shapes. Thus, 
if a known object is within the field of view, precise information about its 
size, location, and orientation can be obtained, for example, from the 
parameters of the least squares ellipse fitted to its boundary points. 

The computational procedures described in this report are totally 
insensitive to image rotation, translation, and scale change. So far as 
is known to the authors, no alternative image processing technique exists 
with a capability of producing extremely accurate object parameter esti- 
mates under such conditions. The algorithms presented are thus felt to 
provide the first feasible method for the generation of very precise 
navigational information from the optical images of known landmarks. 

While all of the results contained in this report relate to two 
dimensional objects, it appears that the basic approach is applicable to 
three dimensional image analysis as well. Specifically, it seems feasible 
that a computational procedure could be developed which would be capable 
of producing a replica of the image of a given three dimensional object 
produced by a particular optical system with a specified spatial relation- 
ship to the object. From such a synthetic image, it ought to be straight- 
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forward to extract boundary points which could then be compared to the 
boundary points of the real scene. Iterative adjustment of the spatial 
parameters, used to generate synthetic images could then be accomplished 
by the nonlinear regression program included in this report so as to op- 
timize the fit of the synthetic image to the real image. Such an exten- 
sion of the present work would permit the use of optically derived 
guidance information in such difficult tasks as automatic orbital rendezvous 
and docking of spacecraft. 

In summary, the ability of a digital computer to extract accurate 
guidance and control signals from optical fields has been established by 
this study. Additional work along the lines indicated by this research 
should eventually produce a very valuable means for a spacecraft or 
robot vehicle to obtain quantitative information regarding its position 
and angular orientation relative to objects within its field of view. 
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APPENDIX I 


GENERATION OF DATA POINTS 


The data points which are presented to the parameter estimation 
algorithm are not physically measured points since no equipment was 
available for this purpose. The entire data acquisition process is 
instead simulated by a digital computer. The following sections briefly 
explain how the data points which lie on the boundary of an ellipse or a 
rectangle are generated, as well as how noisy data points may be gener- 
ated. The subroutine which generated the data points is denoted by DATA. 

Ellipse 


The generation of the data points lying on the boundary of an 
ellipse shall be considered first. An ellipse, as shown in Figure A-l, 
may be expressed analytically by Eqn. A-l. 



Fig. A-l --Ellipse in w, z-reference frame. 


where 


e x w 2 + e 2 z 2 =1 


e x = J— and e 2 = — 

X* 2 x* 2 

W r Z 


(A-l) 


( A - 2) 
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Thus, if one is given the two parameters and e 2 then the corresponding 
ellipse in the w-z plane is completely specified. It is desired to represent 
this ellipse’ s boundary by some finite number of points. For a given 
number of data points, say N, there are infinitely many different ways in 
which these points may be positioned on the ellipse’ s boundary. However, 
it seems quite unrealistic to have the data points very dense on one por- 
tion of the boundary and very sparse, or nonexistent, on the remaining 
portion of the boundary. Perhaps the most realistic situation is for the 
data points to be uniformly distributed on the . boundary of the ellipse. 

This would require the distance between any two adjacent data points to 
be L/N, where L is the length of the boundary of the ellipse. This 
particular distribution of the data points on the boundary of the ellipse 
was not used, however, because of the complex computer programming 
which would be involved and because in practice the physical measuring 
equipment probably would not select the data points in precisely this 
manner anyway. 

The method which was employed to select data points on the 
boundary of the ellipse consists of dividing the w-axis diameter of the 
ellipse into N/2 equal length segments when N data points are desired. 

This, of course, requires N to be an even number. This procedure is 
shown in Figure A-2 for N = 10. 



Fig. A~2--Data points corresponding to ellipse in w, z- 
reference frame. 


The N data points then comprise those points on the boundary of the 
ellipse whose w-coordinates are the same as the w-coordinate s of the 
end points of the segments of the w-axis diameter. 


no 


After the w, z- coordinates of the data points which lie on the 
boundary of the ellipse have been determined, it is necessary to find 
the coordinates of these same data points with respect to the reference 
x, y-coordinate system. These two coordinate systems are shown in 
Figure A-3, 



Fig. A - 3 - - Data points corresponding to ellipse in x,y- 
reference frame. 

Once the x and y-translation (denoted by A and B, respectively) 
of the center of the ellipse and the rotation ( denoted by 0) of the w-axis 
of the ellipse are specified, then the x, y-coordinates of the data point 
having w, z-coordinates (wi,zi) are 

x i = w i cos ^ - z i sin 9 + A (A-3) 

y^ = wi sin 0 + cos 0 + B (A- 4) 

The x, y-coordinate s of the data points are then taken as the 
coordinates of the data points which represent the ellipse whose parameters 
are now to be estimated. 

Rectangle 

The data points which lie on the boundary of a given rectangle 
are generated in a slightly different manner than those of an ellipse. 

Figure A-4 shows a rectangle whose center is at the origin of the w, z- 
coordinate system. 


Ill 




The rectangle 1 s dimension in the w-direction is 2r w , while its 
dimension in the z-direction is 2r z . The dimensions r w and r z may be 
thought of as n radii M of the rectangle. 

The data points are selected such that one quarter of the total 
number of data points lie on each side of the rectangle. This requires 
N to be devisible by four. The data points are further restricted to be 
equally spaced along each side. Therefore, the spacing between data 
points which lie on the vertical boundaries is 8r z /N while the spacing 
between data points which lie on the horizontal boundaries is 8r w /N. 
Once the w, z-coordinates of all the data points lying on the boundary 
of the rectangle are found, their corresponding x, y-coordinate s may be 
determined from Eqns. A-3, 4. 

Noise 


The preceeding discussion has briefly explained how the data 
points corresponding to either an ellipse or a rectangle are generated, 
being given the parameters ej , e 2 , A, B, 9 or r w , r z , A, B, 0. These 
data points fall exactly on the boundary of the appropriate pattern. 

Since there is no error in the coordinates of these data points, they may 
be considered as noiseless data points. 

In a realistic system, however, one would expect that the mea- 
surement points would not exactly overlay the boundary of the pattern 
from which they came. This error may be due to several different 
reasons. For instance, if the field of view has been slightly clouded 
over, or defocused, then the boundary of the pattern is no longer precise 
and the exact coordinates of points lying on the boundary can only be 
estimated. Even if the field of view is clear there is still the possibility 
that the electronic equipment associated with the optical system can com- 
mit errors, be they internal or transmission errors. Furthermore, 
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there is always the quantization error associated with analog to digital 
conversion. All of these errors may be considered as forms of noise. 

Therefore, in order for the artificially generated data points to 
realistically correspond to physically measured data points it is 
necessary to degrade the artificially generated data points by corrupt- 
ing them with some type of noise. A detailed analysis of the physical sys- 
tem would be required in order to know the exact nature of the actual 
noise; i. e. , its distribution and whether it is additive, multiplicative, 
or whatever. In this study no particular physical system was considered; 
therefore, the noise samples were assumed to be additive, statistically 
independent, gaussian noise samples. How well this artificial noise 
resembles the actual noise in a physical system was not considered. 

The gaussian noise was generated by the subroutine GAUSS which 
is in the library of the IBM 360/7 5 at The Ohio State University Computer 
Center. The subroutine permits one to specify both the mean and the 
standard deviation of the gaussian noise samples which it is to generate. 
The subroutine makes use of another library subroutine called EANDU 
which generates uniformly distributed random numbers in the range 0-1. 
The subroutine GAUSS approximates a gaussian random variable by 
adding together twelve uniform random variables, making use of the 
Central Limit Theorem. 

Since it is assumed that the noise is additive gaussian, the numbers 
which are generated by GAUSS are simply added independently to each 
coordinate of the data points. Thus, if the coordinates of a noiseless data 
point are given by (X|,Yj[) , then the coordinates of the corresponding 
noisy data point are ( xi, yi) , where 

x^ = Xi + nj ( A- 5) 

Yi = Y i + n j + 1 ( A-6) 


where nj and nj + \ 


are two consecutive noise samples. 


113 



APPENDIX II 


DISCUSSION OF CRITERION FUNCTIONS 


It is important to recognize the fact that the criterion function 
used in the case where the error function is a linear function of a set of 
parameters ( one step minimization method) is not identical (even to 
within a scale factor) to the criterion function which results when the 
error function is a nonlinear function of a different set of parameters 
( iterative minimization scheme) . 

To be more specific, suppose that the criterion function resulting 
from the linear error function ( Eqn # 3. 53) is denoted by ^ and written 
as 


<f>L.( x, y; p) = e e 
N 

= 2, (p ‘ x i + Pz*iYi + Psyf + P4*i + Ps Yi + 1) 2 
i = l 


(A-7) 


The minimization of is taken with respect to the parameter 
vector p = ( Pi > P 2 » Ps > P 4 > Ps - If p' is that value of p which re suits in 
attaining its unique minimum value, then 


<j>L ( x, y; p" ) = min <t> L ( x, y; p) 


( A-8) 


The ellipse parameter vector c = ( e x , e 2 , A, B, 0) is then found 
from Eqn. ( 3. 42) , that is. 


c 


» ftp*) 


( A-9) 


The criterion function resulting from the nonlinear error function, 
denoted by may be written as 


4> N (x^ yi?Cc) ) 


-TV 
e e 



p ix2 + p 2 xiyi + P sYi + P 4 x i + P 5 Yi + P6 + 1) 2 

( A- 10) 
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where the vector p = ( p 1 » P 2> Ps» P 4 * P 5 > P6) is a nonlinear function 
of the ellipse parameter vector c? with which the minimization of ^ is 
taken. The vector ]?(<?) is given by Eqns. ( 3. 10 - 3. 15). 


Note that 4 ^ may be factored 
N 


4^ = ( P6 “ 1 ) [ ■ ft * 1 x i + -^- 2 — x iYi + — P ’ 3 ■ y\ + P * x i + - 2-5 . yi + 1 ] 

. ^ P 6 “ 1 P6“l P6“l P6 “ l P6-1 

1 = 1 


= ( P6 - 1 ) 2 4>L 


(A- 11 ) 


Since p* is a function of p* ( Eqn. 3. 39) » the factor ( P6 _1 ) 2 may 
be represented by some function, g( p) , to give 


4> N ( x >y;p) = g(p) 4> T .( x »y5p) 


( A - 1 2) 


If p*^ is that value of p for which 4 >n attains its minimum value. 


denoted by then 

❖ 


<I> N ( X, y ; p" ‘ ) = min <(> N ( x, y; p) 


( A- 1 3) 


and so 


4> N ( x, y; p 


) 1 g(p ) 



( A -1 4 ) 


One should note that if 4 >n is divided by the factor ( p& -1 ) 2 before 
the minimization with respect to c* is taken, then the two criterion func- 
tions would be identical and both minimization techniques would yield the 
same value for the minimizing parameter vector ( assuming no boundaries 
are encountered). 


5jC jfc 

A comparison of (j^ and is made in Table A. The right hand 
side of Eqn. ( A-14) is also tabulated, being denoted by cjT. The values 
of the criterion functions are those that were obtained by using the two 
minimization schemes on an ellipse which is characterized by 


c 0 = 


' 0. 25 ' 
1 . 00 
1 . 00 
- 2 . 00 
_0. 50. 


( A- 1 5) 


The parameter vector estimates which correspond to these values for the 
criterion functions are shown in Tables 1 and 2 in Chapter IV. 
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It is interesting to note that in every case 
as Eqn. (A-14) implies. 


* ^ 


in Table A-I, 


0 . 0 
0. 0 
0 . 0 
0. 0 

0. 1 
0. 1 
0. 1 
0. 1 

0 . 2 
0 . 2 
0 . 2 
0 . 2 

0. 3 
0. 3 
0. 3 
0. 3 

0. 4 
0. 4 
0. 4 
0. 4 

0. 5 
0. 5 
0. 5 
0. 5 


TABLE A-I: COMPARISON OF CRITERION FUNCTIONS 


^L 


<l> 


4 > 


=s= 

N 


10 

0. 1564 x 10“ 7 

0. 2496 x 10“ 6 

20 

0. 2379 x 10“ 7 

0. 3795 x 10 -6 

50 

0. 6631 x 10“° 

0. 1060 x 10“ 4 

100 

0. 8263 x 10“ 6 

0. 1319 x 10“ 4 


0. 6555 x 10” 11 
0. 1846 x 10' 10 
0. 7957 x 10" 10 
0. 4480 x 10“ 10 


10 

0. 1757 x 

10 

20 

0. 1747 x 

10 

50 

0. 8430 x 

10 

100 

0. 1697 x 

10 


0. 3219 x 10° 
0. 3099 x 10 
0. 1313 x 10 1 
0. 2710 x 10 1 


0. 2911 x 10° 
0. 2913 x 10° 
0. 1191 x 10 1 
0. 2449 x 10 1 


10 

20 

50 

100 


0. 2264 x 10” 1 
0. 1611 x 10° 
0. 2627 x 10° 
0. 5216 x 10° 


0. 3992 x 10° 
0. 2566 x 10 1 
0. 5528 x 10 1 
0. 1164 x 10 2 


0. 3467 x 10° 
0. 1687 x 10 1 
0. 3989 x 10 1 
0. 8214 x 10 1 


10 

20 

50 

100 


0. 9620 X 10' 1 
0. 3844 x 10° 
0. 7211 x 10° 
0. 1881 x 10 1 


0. 2193 x 10 1 
0. 6326 x 10 1 
0. 1151 x 10 2 
0. 3245 x 10 2 


0. 1232 x 10 1 
0. 2912 x 10 1 
0. 6331 x 10 1 
0. 1480 x 10 2 


10 

0. 2027 x 

10° 

0. 9334 X 

10 

20 

0. 2268 x 

10° 

0. 1817 x 

10 

50 

0. 1097 x 

10 1 

0. 3412 x 

10 

100 

0. 1 544 x 

10 1 

0. 1190 x 

10 


0. 5916 x 10° 
0. 6042 x 10 1 
0. 9126 x 10 1 
0. 1520 x 10 2 


10 

0. 6232 x 

10 

20 

0. 4839 x 

10 

50 

0. 1257 x 

10 

100 

0. 3031 x 

10 


0. 7164 x 10° 
0. 1388 x 10 2 
0. 7153 x 10 2 
0. 1048 x 10 3 


0. 5555 x 10° 
0. 4421 x 10 1 
0. 8891 x 10 1 
0. 2317 x 10 2 
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APPENDIX III 


COMPUTER PROGRAM 


A listing of the Fortran program which was used to estimate the 
five parameters associated with an ellipse is given in this appendix. The 
program has been broken down into a main program along with several 
subroutines, each of which has a specific function. Each subroutine is 
briefly discussed in the following paragraphs. 

MAIN Program 

The main program performs three functions. The first function 
is to read all the required input information for the overall program. 
Secondly, the main program calls the various subroutines in the correct 
sequence such that the iterations for the estimates of the parameters are 
correctly performed. Finally, the main program writes out the input 
information as well as the best estimate for the parameter vector. 

The main program that is listed in this appendix is the one which is 
used in the estimation of the parameters of an ellipse. The main program 
which is used for the estimation of the parameters of a rectangle is iden- 
tical to the listed main program except that one statement is added which 
relates the size of the fitted ellipse to the size of the corresponding rec- 
tangle, This scale factor is discussed in Chapter V. 


The main program requires the following inputs: 


NPAR the number of parameters which are to be 

estimated. 


NPOINT the number of data points which lie on 
the boundary of the unknown pattern. 

NTRIAL the total number of independent local 

minimizations of the criterion function 
with respect to the parameter vector. 
NTRIAL > 1. 
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MSMSQ 

MR AND 

NSET 

CE 

MINPAR 

MAXPAR 

SIGMA 

AVE 

FSTOP 

EPHI 

EC 

EGRAD 


the maximum number of times which the 
LOCMIN subroutine may call the SUMSQR 
subroutine, 

the number of independent parameter vectors 
which are randomly selected by the RANSER 
subroutine. 

the number of patterns whose parameters 
are to be estimated. 

the initial guess for the unknown 
parameter vector. 

the vector corresponding to the lower 
limits for the parameters. 

the vector corresponding to the upper 
limits for the parameters. 

the standard deviation of the gaussian 
noise which is added to the simulated 
data points generated in the DATA 
subroutine. 

the average value of the gaussian 
noise associated with SIGMA. 

the upper limit for the absolute value 
of F. The program is terminated if If I 
becomes larger than FSTOP. 

the lower limit for DPHI, The program 
is terminated if DPHI becomes smaller 

than EPHI, where d± = ^ ~^( c i + 1 1 

* « Ci) 

the lower limit for DC. The program is 
terminated if DC becomes smaller than 
EC, where 

d c = “ 

c i c i 


the lower limit for the squared magnitude 
of the gradient vector. The program is 
terminated if |vcK Ci )| becomes smaller 
than EGRAD. 
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EBDBY 


a constant which is used in the GRPEEX 
subroutine to prevent division by zero. 


DATA Subroutine 

The purpose of the DATA subroutine is to artificially generate the 
data points which lie on the boundary of either an ellipse or a rectangle. 
Although the DATA subroutines corresponding to both an ellipse and a 
rectangle are shown in the listing, only one of them is included in the 
program when it is actually used. Appendix I gives more details as to how 
the data points are generated and how simulated noise is added to them. 

SUMSQB Subroutine 


The SUMSQB subroutine simply evaluates the criterion function 
for a specific value of the parameter vector. It also has an instability 
indicator, KX, which is set to one if the absolute value of F exceeds 
FSTOP. 

LOCMIN Subroutine 


The LOCMIN subroutine performs a local minimization of the 
criterion function with respect to the parameters. It does this by calling 
the next three subroutines. It also checks the various criteria for ter- 
minating the program. 

BEGBES Subroutine 


The BEGBES subroutine evaluates the criterion function ( PHI) , 
the gradient of the criterion function ( GBADP) , and the Gauss-Newton 
parameter change vector ( BETA) for a specified parameter vector ( C) 
which is supplied by the LOCMIN subroutine. A library subroutine ( MINV) 
is used for matrix inversion. 

GBASEB Subroutine 


The GBASEB subroutine is called only when the Newton-B aphson 
method is used to determine the next value for the parameter vector. The 
GBASEB subroutine finds the optimum binary scale factor by which to 
multiply Acjp 
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G RPR EX Subroutine 


The GRPREX subroutine is called only when the full Newton-Raphson 
. _ tfc-) V^ci) 

step (Ac.: = - — ) violates a range constraint. It then projects 

lw ?i) 1 2 

the gradient onto the constraint surface, after which the GRASER subrou- 
tine finds the optimum binary scale factor for this projected gradient. 

The GRPREX subroutine has an output variable, KEXIT, which when set 
to one indicates that the parameter vector is on a constraint boundary of 
the constraint region, R. 

RANSER Subroutine 


The RANSER subroutine selects a given number ( MRAND) of 
parameter vectors randomly, using a uniform distribution, and determines 
that parameter vector which yields the smallest value for the criterion 
function. This parameter vector is then used as the initial guess for 
another local minimization. The RANSER subroutine uses a library sub- 
routine, RANDU, for its uniform number generator. 

A listing of the two library subroutines, RANDU and GAUSS, which 
generate numbers possessing uniform and normal distributions, respec- 
tively, is given at the end of the program listing for the sake of complete- 
ness. 
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SUBROUTINE RANDU ( IX, IY, YFL) 

IY = IX * 65539 

IF (IY) 5, 6, 6 

IY = IY + 2147483647 + 1 

YFL = IY 

YFL = YFL * . 4656613 E-09 

RETURN 

END 


SUBROUTINE GAUSS (IX, S, AM, V) 

A = 0. 0 

DO 50 I = 1,12 

CALL RANDU (IX, IY, Y) 

IX = IY 
A = A + Y 

V = ( A- 6. 0) * S + AM . 

RETURN 

END 
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DIMENSION COUO),CE( iO»tCll 10) 

REAL M l NPAR ( 10),MAXPAR( 10) 

common/comi/fstcp 

CUMMCN/CUM3/NGAUS 
NG AU S = 317578125 
K*l 

1 CALL SCL0K1 

READ 15,5001) NPAR.NPOINTrNTRlALrMSMSOrWRANDfNSET 

5001 FORMAT (615) 

READ (5,5002) C CO I 1 ) , I ■ 1 ,NPAR ) , ICE! I) , I - 1 , NP AR ) , < M I NP AR ( I ), I-l.NPA 
IK I , (MAXPAR ( I ), 1 = 1, NP AH), SIGMA, AVE, ESTOP, FPH I , EC ,E GRAD , EBDR Y 

5002 FORMAT (5F14.6) 

WRITE (6,5003) 

5003 FORMAT (•!• /////////////////• 

1 NONLINEAR SYSTEM PARAMETER ESTIMATION*//////) 

WRITE (6,5004) 

5004 FORMAT (• NP AR NPOINT NTRIAL 

1 MSMSQ MR AND NSET*/) 

WRITE (6,5005) NPAR , NPO I N T , NT R I AL , M SM SQ , MRANO , NSE T 

5005 FORMAT I I 1 0 , 5 ( 1 OX , I 10 ) ) 

WRITE (6,5006) 

5 C06 FORMAT (///* SIGMA AVE FSTOP 

L EPHI EC EGRAD EBDRY • / ) 

WRITE (6,500 7) S I GMA , AVE , F S TOP , E PH I , EC , EGR AD , E BOR Y 
500 7 FORMAT { T 1 2 * F6 . 3, T 28 , F 5 . ?• T 38, F 1 3. 1 , T 55 , F 1 2 . 1 0 , T69 ,F i 2 . 10 . T8 3 , F 1 2 . 
110,1 100, F 12. 10) 

WRITE (6,5008) 

5008 FORMAT (///* I C0( I ) 

1 CE ( I ) MINPARU) MAXPAR ( I ) * / ) 

WRITE (6,5009) ( I , CO ( I ) , C E ( I ) ,M INPAR ( I ) .MAXPAR ( I ) , 1*1, NPAR ) 

5009 FORMAT t I 20 , T 36 , F 8 . 4 , T 56, F 8 . 4 , T 76 , F 8 . 4, T 96 , F 8 . 4 ) 

NSMSC=0 

CALL DATA ( CO , S I GMA , AV E , NPO I NT ) 

CALL SUMSQR < CE , PH I , K X , NSMSQ . NPO I NT , NP AR ) 
f F ( K X ) 3,3,2 

2 WRITE (6,5010) 

5010 FORMAT (//////////• THE INITIAL PARAMETER ESTIMATES PRODUCE AN UNS 
[TABLE RESPONSE. DESCENT TO A MINIMUM WILL NOT BE CARRIED OUT.*) 

GO TC 4 

3 CALL LCCMIN ( CE , PH I , M I NPAR , MAXP AR , EPH I , EC , EGR AD , EBDR Y , MSMSQ , NPO I N T 
l.NPArt) 

IFINTRIAL — 1 ) 8,8,4 

4 IY=123456789 

DU 7 1=2, NTRIAL 

CALL RANSFR ( C I , M I NPAR , MA XP AR , MRAND , I Y, NPO I NT , NPAR ) 

CALL LUCMIN ( C 1 , PH I LOC , M I NPAR , M AXP AR , E PH I , EC , EGRAD , E BOR Y , M SMSQ , NPO 
1INT.NPAR) 

I F ( PH I -PH I LOC ) 7,7,5 

5 PH I = PH I LOC 

DO 6 J = 1 , NPAR 

6 CE ( J I *C 1 ( J ) 

7 CUNTINUE 

8 WRITE (6,5011) 

5011 FORMAT (*1* ////////////////• 

1 FINAL PARAMETER E ST I M A TE S * / / / / / * 

2 C( 1 ) C ( 2 ) C ( 3 ) CK) 

3 C ( 5 ) • / ) 

WRITE (6,5012) ( CE ( I ) , I * 1 , NP AR ) 

5012 FORMAT ( T 1 5 , F 1 3 . 7 , T 35 , F 1 3 . 7, T 55 , F l 3 . 1 , T 75 , F 1 3 . 7 , T95 , F l 3 . 7 ) 

WRITE (6,5013) PHI 

5013 FORMAT ( ////////////////* 

l SUM-SQUARED ERROR * *,E15.8) 

T I M E = RC L OK 1 ( 1.) 
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THE TIME REQUIR 


WRITE I6,50 1 *) TIME 
501* FORK AT (//////////• 

IED TC EXECUTE THIS PROGRAM WAS*,F7.2,* SECONDS. •) 

9 IFCK-NSET) 10,11,11 

10 K*K4l 
GO TC l 

11 CONTINUE 
STOP 
END 


THIS SUBROUTINE IS USED WHEN DATA POINTS LYING ON THE BOUNDARY OF 
AN ELLIPSE ARE REQUIRED. 

SUBROUTINE DAT A I C , S IGMA , AVE • NPO INT ) 

DIMENSION XOI 1001* YO I 1001, XDC 100) tYDUOO 1,01 101 
CQMMCN/C0M2/ X t 1 00 ) , Y I 100 ) 

C0MM0N/C0M3/NGAUS 
DATA IZ*/* •/ 

WRITE 16,5001) 

5001 FORMAT (MENTER DATA SUBROUTINE.*//) 

ei*cm 

E 2*C ( 2 ) 

A*C ( 3 ) 

6=C(*) 

TH = C I 5 ) 

RX-SCRTM ./El ) 

XINI=?.*RX 
XUI 1 ) *— RX 
YOI l ) *0 . 

1=2 

1 XOt I )=-RX*I*XlNT/NPCINl 
XU( l*l)=XO( I ) 

YOI I ) = SQR T 1 I 1 • /E2 ) * I 1.-E1 *X0 I I ) **2 > ) 

YOI l ♦ 1 ) *— YO I I ) 

1*1*2 

IF U-NPOINT 1,2,2 

2 XUlll-RX 
YOI I )=0. 

WRITE 16,5002) CIl),CI2) 

5C02 FORMAT UX,* THE FOLLOWING POINTS ARE THE COORDINATES OF AN ELLIPS 
IE HAVING NO DISPLACEMENT FROM THE ORIGIN*/* AND NO ROTATION ABOUT 
2 THE ORIGIN. THE EQUATION OF THF ELLIPSE IS GIVEN BY Ei*X*X ♦ E2*Y 
3* Y * 1*/* WHERE El =*,F10.5,* AND E2 **,F10.5//) 

N S* NPO I NT— 1 

WRITE 16,500 3) I I Z ft , XOI I ) , YOI I ), 1 = 1, NS), IZ ft , XOI NPO INT ) , YO I NPOI NT ) 
6003 FURMAT 1 * I A 1 , • I • , E 1 2 . 5 , • , • , E 12 . 5, • ) • ) ) 

DO 3 I = 1 , NPO INT 

XDC I )*XOt I >*CGS I TH)-YCJ< I ) *S IN 1 1 H ) *A 
YDI I ) *X0( l )*SlN( TH)*YOI I ) *COS I TH ) *B 

3 CONTINUE 

DEGREE * I 1A0./3. 1 A 1 59265 36 ) * TH 
WRITE I 6 , 5004 ) C(3),C(4), DEGREE 

500* FORMAT I//2X, *THE FOLLOWING POINTS ARE THE COORDINATES OF THE PREVI 
20US ELLIPSE WHICH HAS NOW BEEN DISPLACED A =*,F10.5/* UNITS IN THE 
3 X-DIRECI 1 ON , B =*,F10.5,* UNITS IN THE Y-D1RECT10N, AND ROTATED b 
3Y THETA =*,F10.5,* DEGREES.*//) 

WRITE (6,600 5) I I Z ft , XD I I ) , YD (I) • I * 1 , Nft ) , l Z ft , XD I NPO I NT ) , YD I NPO I NT ) 
5C05 FORMAT I * I A l , • I • , E 1 2 . 5, • • • , E 12. 5, • ) • ) ) 

DO * 1=1, NPO INT 

CALL GAUSS I NGAUS , S I GMA , AVE , V ) 

XI I ) =XD I I)*V 

CALL GAUSS I NG AUS , S I GMA , AVE , V ) 

* Y I I ) = YO I I ) ♦ V 


123 



o n 


5006 

5C07 

5008 

5C0 1 


1 

2 

3 

4 

5 

6 
7 


Write «6»5006) ave*sigpa 

FGRPATI//2X* 'THE following points are the previous points with gau 

1SSIAN NOISE ADDED. THE MEAN OF THE NOISE IS c tEl2.5/ a AND THE STAN 
2DAKC DEVIATION OF THE NOISE IS **E12.5//) 

WRITE 16 1 50071 MZ$»X< I ) * Y (II • I* 1 tNS ) * I Z S* X I NPOI NT ) # Y I NPOI NT ) 
FORMAT ( 4» ( Ai * ' C , tEI2«5t # * # tE12«5t f ) • ) ) 

WRITE (6 • 5008 ) 

FORMAT l/' EXIT CATA SUBROUTINE 1 ) 

RETURN 

END 


THIS SUBROUTINE IS USED WHEN DATA POINTS LYING ON THE BOUNDARY OF 
A RECTANGLE ARE REQUIRED. 

SUBROUTINE DAT A ( C t S IGM A , AVE * NPO IN T ) 

DIMENSION CIS) ,XOUOO)*YOI 100)*XDI 100) * YDC 100)* XII 100) ,Y1 l 100) *X2( 
1 100) * Y2U00) ,X3U00) , Y3( 100) ,X4( 100) * Y4| 100) 

CUP MON/ CUP 2/ X ( 100) * YI 100) 

COMMON/ CUM 3/ NGAUS 
DATA 11%/* */ 

WRITE 16,5001) 

FORMAT I’lENTER LATA SUBROUTINE.*//) 

XMAX=1.0/ SQRT ( C I 1 ) ) 

YMAX = 1.0/S0RT(CI2) ) 

A*C I 3 ) 

B = C ( A ) 

TH.CI5) 

C TH*CUS ( T H ) 

STH*S IN ( TH) 

Y M | N=— YMAX 
XMIN^-XMAX 
NP*NPOINT/A 

X INT=XMAX-XMIN 
YINT =YMAX— YMIN 
\=NF*M 
DU 1 I * 1 , N 
XIII ) = X MAX 

Ylli )=YMIn*< 1-1 )*Y INT/NP 

DO 2 I = 2 , NP 

X2( I )*XMAX-( I- l )*X INT/NP 
Y 2 C I ) * YMAX 
DU 3 1=1* N 
X 3 I I ) = XM IN 

Y 3 I I ) * YMAX— I I- 1 )»Y INT/NP 
DO A 1*2, NP 

X4m=XPtN+( 1-1 )*X INT/NP 
Y<* ( I ) *YM I N 
DO 5 I =1 ,N 
X0( I )*X1 1 I ) 

YO ( I ) *Y l II ) 

N*NP* 2 
NN* 2 *NP 
DO 6 I *N, NN 
X0( I ) =X2 I I-NP ) 

YOU ) = Y2 ( I-NP) 

N = 2 *NP* l 
NN*3*NP* l 
DO f I=N,NN 
XO ( I ) *X 3 I I — 2 *NP ) 

YO ( I ) * Y 3 ( I — 2*NP ) 

N*3*NP*2 

nn=a*np 

DU 8 I =N» NN 
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X0| 1 )rX A( i-3*NP> 

YOU )=YM I-3*NP) 

WRITE 1 6 , SCO 2 ) XINT.YJNT 

5002 FORMAT (2X,*THr FGLLGWlNG POINTS LIE ON A RECTANGLE HAVING NO OISP 
1LACEMENT FROM THC ORIGIN ANU NO ROTATION ABOUT THE ORIGIN. •/' THE 
2REC F ANGLE IS CENTERED ABOUT THE ORIGIN AND IS*,F6.2,* UNITS IN LE 
BNGTH IN THE X-C IKECT ION AND * , F h. 2 , * UNITS*/* IN LENGTH IN THE Y-D 
4IRECFICN. •) 

Ni~N PO I NT- I 

WRITE (6,5003) I I Z S • XO I I ) , YO ( I ) , l - 1 , NS ) , I 2 S , XOI NPO INT ) ,Y0 INPOINT ) 

5C03 FORMAT K(Ai,M*,E12.5,'t*,El2.5,*)*n 
DO 9 1*1 * NPO INT 
XL) ( I I =X0 ( I )<=CTH-YO( I J « S T H t A 

9 Y0( I )=X0( D + STH + YOI 1)*CTH+B 
WRITE I 6 « 5004 ) A , 8 » T H 

5004 FORMAT I//2X,*THE FOLLOWING POINTS ARE THE POINTS ON THE PREVIOUS 
iRECTANGLE WHICH HAS NOW BEEN T R ANSL A T EO ' , F6 . 2 , * UNITS IN THE*/' X 
2-DIRECT ION, • ,F6.2, • UNITS IN THE Y-OIRECTION, AND ROTATED BY * , F6 • 
32,* RADIANS.*) 

WRITF (6,500 5) ( I £ t • XOCfl • YD ( I) • T= 1 • Nt ) « I Z $ , XD ( NPO I N T ) , YD (NPO I NT ) 

5C05 FORMAT (A ( A1 , • I • , E 1 2 . 5 , • • * • E 1 2 . 5 , * ) * ) ) 

DO 10 1*1, NPC (NT 

CALL GAUSSINGAUS, SIGMA, AVF,V) 

xm=xcm*v 

CALL GAUSSINGAUS, SIGMA, AvE.V) 

10 Y I I ) =YD ( I ) ♦ V 

WRITE (6,5006) AVE, SIGMA 

5 C06 FORMAT I//2X, *THE FULLOWlNG POINTS ARE THE PREVIOUS POINTS WITH GAU 
1SSIAN NOISE ACDEO. THE MEAN OF THE NOISE IS *,EI2.5/* AND THE STAN 
2DARU DEV I A T I UN CF THE NOISE IS *,E12.5//) 

WRITE 16,6007) (IZS,X(I),Y(I),I=1,NS), I Z l , X ( NPO I N T ) , Y I NPO I N T ) 

5C07 FORMAT IMAl,*t*,E12.5,*,*.E12.5,*)*)) 

WRITE I 6 , 500 8 ) 

5008 FORMAT (/• EXIT DATA SUBROUTINE*) 

RETURN 

END 


SUBROUTINE SUMSCR ( C.PHl ,KX,NSMSQ, NPOiNT ,NPAR ) 

DIMENSIUN C( 10) ,FI 100) 

COMMCN/COM1/FSTCP 
C0MMCN/CUM2/XI ICO) , Y( IOC) 

NSMSC=NSMSC + 1 
KX = 0 
PH I =0 • 

E 1 =C ( 1 ) 

E 2 = C ( 2 ) 

A- C ( 3 ) 

B=C (A) 

TH=C ( 5 ) 

C T H = CUS I f H) 

STh=SINI TH) 

DU 2 I = l • N PC I NT 

F( I ) *£ 1 * ( (XI I)-A)*CTH*(YI I )-B)+STH)**2*E2*l-(X( I)-A) *STH* ( Y ( I ) -B ) * 
1C TH ) ♦ *2— 1 ■ 

IF ( ABS ( F ( I ) )-FSTOP ) 2,1,1 
1 KX=I 

WRITE (6,500?) ( C ( K ) , K= 1 , N.P AR ) 

500? FORMAT I* THE SYSTEM IS UNSTABLE FOR THESE PARAMETER VALUES. THE S 
1UM-SCUARED ERROR WILL NOT BE EVALUATED. THE PARAMETERS ARE * / * C C 1 ) 
2 =*.FI5.8,* C ( 2 ) =■ , E 15. 8, • C 1-3 I =*#El5.fl,* C l 4 ) «*,E15.8,* C(5) * 
3*,E15.R) 

GO TO 3 
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2 PHI « pH I ♦ F( i )*F( I) 

l CONTINUE 

RETURN 
END 


SUBROUTINE LOCK IN ( C I , PH I , H I NPAR , MAXP AR • EPH I , EC • EGRAD , EBDR Y , HSMSQ , 
INPOINT, NPAR) 

DIMENSION U 10) ,CI ( 10 ) , 6RA0P ( 10) .BETA I 1 0 ) , DEL TAC ( 1 0 ) 

REAL MINPARt 10) ,MAXPAR( 10) 

NSMSC*0 
NS = 0 
ND = 0 

WRITE (6,5001) 

5001 FOR HAT ('IENTER LOCMIN SUBROUTINE.*) 

1 CALL REGRES ( C 1 , GR ACP , B E T A , PH I , NPO I N T , NPAR ) 

PHI I N T = PH I 

csc=o 

DU 2 1=1, NPAR 

2 csq=csq+ci ( i >*cim 

GRACPS=0. 

DO 3 1=1, NPAR 

3 GKACPS = GRADPS + GRADP( I ) * GR AOP ( I ) 

IF (GRAOPS-EGRAD ) A , 5 , 5 

A WRITE (6,5002) GRADPS 

*5002 FORMAT (/• THE GRADIENT CONDITION IS SATISFIED. THE GRADIENT MAGNI 
ITUDE S0UAR6D IS GRADS = • , E 1 5. 8 ) 

GO TC 32 

5 bU 6 1=1, NPAR 

6 C ( I ) =C1 ( I > ♦BET A ( I ) 

DO ti 1=1, NPAR 

IF (C ( I )-M. INPAR ( I ) ) 20,20, 7 

7 IF tC(!)-PAXPAR( II) 8,20,20 

8 CONTINUE 

CALL SUMSCJR (C,PhIREQ,KX,NSMSO,NPOINT,NPAR) 

IF(KX) 9,9,21 

9 I F ( PH I REG-PH I ) 10,22,22 

10 00 11 1=1, NPAR 

u c i »=cm 

11 OF L T AC ( I ) = BE l A ( I ) 

WRITE (6,5003) PH I REG 

5C03 FORMAT (/• THE FULL GAUSS-NEWTON STEP YIELDS A SMALLER VALUE FOR 

1 PH I WITHOUT VIOLATING*/* ANY CONSTRAINTS. THE NEW VALUE FOR PHI IS 

2 PHI = *,E15.ei 
WRITE (6,5004) 

500* FURMAT (/• I C1U)*/) 

WRITE (6,5005) ( I , C 1 I I ) , I - 1 » NPAR ) 

5005 FORMAT U10.E20.8) 

UPHI = (PHI-PHIKEG )/PHl 
PHI =PHI«EG 
GU TL 13 

12 UPHI =(PHI-PHIGRA)/PHI 
PH I = PH I GR A 

13 IFIDPHI-EPHII 14*14,15 

14 WRITF (6,5006) CPH1 

5C06 FURMAT (/• THE NORMALIZED SUM-SQUARED ERROR CHANGE CRITERION IS SA 
ITISFIED. UPHI = *,E15.8) 

GO TU 32 

15 OELCSQ-O. 

DO 16 1*1, NPAR 

16 DELCSQ = DELCSC*DELTAC( I ) *DEL T AC l I ) 

DC*OFLCSQ/CSQ 

IF (DC-EC) 17,17,18 
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17 WRITE (6,5007) CC 

5C07 FORMAT (/• THE PARAMETER CHANGE ERROR CRITERION IS SATISFIED. DC * 
1 • , E 1 5 • 8 ) 

GO TC 3 2 

18 IF(NSMSQ-MSMSQ) 1,19,19 

19 WRITE (6,5008) 

5008 FORMAT (/• THE LOOP COUNT CRITERION IS SATISFIED.*) 

GO TO 32 

20 WRITE (6,5009) 

5C09 FORMAT (/• THE GAUSS-NEWTON PARAMETER CHANGE VECTOR VIOLATES THE 
1 RANGE CONSTRAINTS. •/' A NEWTON-R APHSON STEP WILL BE TRIED.*) 

GO TC 23 

21 WRITE (6,5010) 

5C10 FORMAT (/• THE GAUSS-NE WTCN STEP YIELDED AN UNSTABLE VALUE FOR PH 
II. THEREFURE, •/• A NEWTON-R APHSON STEP WILL BE TRIED.*) 

GO TC 23 

22 WRITE (6,5011) 

5011 FORMAT (/• PHI OBTAINED FROM THE GAUSS-NEWTON STEP WAS GREATER TH 
IAN PHI BEFORE THt STEP WAS TAKEN. ■/• THEREFORE, A NE WTON-R APHSON S 
2TEP WILL BE TRIED.* ) 

23 DO 24 1=1, NPAR 

DE L T AC ( I ) =- ( PH 1 *GR AC P ( I ) >/GRADPS 
2* Cl I )=C1 ( I ) ♦CELTACI I ) 

WRITE (6,5012) 

5012 FORMAT (• I HE FULL NEW T CN-R APHSON STEP IS*//* 1 DE 

ILTACII) CII)*/) 

WRITE (6,5013) ( I, DEL T AC ( I ),C( I ), 1= l,NPAR) 

5013 FORMAT I I 1 0 , 5X , E I 5 . 8 , 5X , E 1 5. B ) 

DO 26 1=1, NPAR 

IF (C( I )— M I NPAR ( I ) ) 30,2 5,25 

25 I F I C I I )-MAXPAR< 1 ) 1 26,26,30 

26 CUNT1NUE 

27 CONMNUF 

B 1 NSC L = 2 • * *NS 
UU 28 1=1 , NPAR 

28 DELTACf I )=DELTAC( I l/HINSCL 

CALL GRASER ( PH I I NT , C 1 • NS , DEL T AC , K X , NSM SQ , NPO I N T , NPAR ) 

phigra=phi int 

NT=NS 

NS*NC 

ND = NT 

IF (KX) 29,29,32 

29 WRITE (6,5014) ND , NS 

5014 FORMAT (/• T HF LPT I MUM DELAY SCALE FACTOR FOUND FRUM THE GRASER S 
1UHR0UTINE WAS N = *,I2,* AND THE INITIAL DELAY SCALE FACTOR TO BE* 
2/* USED FUR THE NEXT GRASER SUBROUTINE (IF IT IS CALLED AGAIN DUR I 
)NG THIS LOCMIN SUBROUTINE) IS N = *,I2) 

WRITE (6,5015) PhIGRA 

5015 FORMAT (/• THE GRASER SUBRUUTINF FOUND PHI = *,E15.8,* AND THE PA 

1RAMETERS TU BE*//* I C1(I)'/) 

WRITE (6,5016) (I,C1(I),I=1,NPAR) 

5016 FORMAT ( I 10 , 5X , E20 • 8 ) 

GO TC 12 

30 WRITE (6,5017) 

5017 FORMAT (/• THE NE W T ON- R APHSON STEP VIOLATES THE RANGE CONSTRAINTS 
1.'/* GRADIENT PROJECTION AND EXTRAPOLATION FOLLOW.*) 

CALL GRPREX ( C l , DE L T AC , M I NP AR , MAXP AR , GR AOP , KE X I T , E BDR Y , NPAR ) 
IMKFXIT) 2 7,27, 31 
>1 WRITE (6,5018) PHI 

5 1 1 0 FORMAT (/• THE CUNSTRA I NED MINIMUM IS PHI * *,E15.8,* AND THE PARA 
1METERS ARE*//* I C 1 I I) */ ) 

WRITE (6,5019) ( I, CM I ), 1*1, NPAR) 

501? FORMAT f I 1 0 , 5X , F 1 5 . 8 ) 
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32 CUNTINuE 

WRITE (6.5U20) NSMSQ 

5020 FORMAT!//* EXIT LOCMIN SUBROUTINE. NSMSQ « "«I*> 
RETURN 
END 


SUBROUTINE RFGRES ( C t GRACP, BET A, PH I • NPO I NT, NPAR ) 

DIMENSION C ( 1 0 ) » GRADP ( 10 ) * BET A ( 10 ) » F 1 100) « l ( 100, S j , S ( 5,5) ,L< 5) ,M< 
15) 

C0MHCN/C0M2/XI IC0),Y( 100) 

WRITE (6.5001) 

5001 FORMAT (/• ENTER REGRES SUBROUTINE.*) 

E l =C ( 1 ) 

E 2 *C ( 2 ) 

A = C< 3) 
d*C(4) 

TH»C ( 5 ) 

CTH-COSl TH) 

STH=SIN(TH) 

00 l 1 = 1. NPO I NT 

/(I,n*imi)-A)*CTHMYII)-B)*STH)**2 
/I I . 2 ) = ( - I X I I ) — A)*STH+(Y( I )-B)*CTH)**2 

Z< I »3)=2.*(E1*CTH**2*E2*STH**2)*(— X( I )+A) ♦ 2.*(El-E2)*STH*CTH*(-Y 

l ( I 1 ♦ tl I 

1 ( I ,4 )=2.*( 1 1-F2)*STH*CTH*|-X( I )*A) «■ ? . * ( E 1 * STH**2^E 2 *C TH**2 ) * ( - Y 

i m*B> 

ZM,5)=2.*(E1-E2)*(l-(X(I)-A)**2*(Y(I )-B)**2)*CTH*STH*( X( I ) - A ) * I Y ( 
II )-B) + (CTH**?-STH**2) ) 

1 F( I )=£!♦( ( X ( I )-A)*CTH+(Y( I ) - B ) * STH ) * *2 ♦ E2* ( - ( X I I ) - A ) *STH* ( Y (I ) -B ) * 


ICTH ) # *2— 1 . 
WRITE (6,5002) 
5002 FURMAT (//• 


Z f I , 1) 

Z ( 1,2) 

Z 

1(1,3) 

Zf 1,41 

ZC 1,51 

F ( I ) * / ) 



WRITE (6.500)) ( I , Z ( I , 1 ) , Z ( I , 2 ) • Z ( I,3)«Z(I,4),Z( 1.5). F(I). 1 = 1. 

1NPT2 ) 

5 C03 FORMAT (I 5 , 5 X , 6E 1 8 . 8 ) 

DO ? 1=1 .NPAR 
DU 2 J=1 .NPAR 
S( I . J ) =0 • 

DO 2 K = 1 . NPC I NT 

2 S( I » J ) =S ( I , J I + Z ( k, I ) *Z ( K , J ) 

WRITE ( 6 . 5004 ) 

5004 FORMAT (//* I S(I.l) S(1.2) S 

HI.)) S( 1,4) S( 1.5)*/) 

WRITE (6.5005) (I.S(I.1).S(I.2).S(I V 3)»S(I.4) V S( I » 5 ) » I = 1 .NPAR ) 

5005 FORMAT ( I 5 . 5 X . 5 E 1 8 . 8 ) 

CALL MINV(S.NPAR.D.L.M) 

WRITE (6.5006) 

5006 FORMAT (//• THE FOLLOWING ARRAY CONTAINS THE ELEMENTS OF S INVERSE 


I. HOWEVER, THEY 

WILL 

ST ILL 

BE DENOTED BY S.'//* 

I 

2(1,1) 

3 S f I * 5 J f / ) 

S( I 

t 2) 

S( I « 3) 

S( 1,4) 


WRITE (6,500 5) ( I , S ( I , 1 ) . S ( l , 2 ) , S ( I , 3 ) , S ( 1 . 4 ) , S ( I , 5 ) , I « 1 , NPAR ) 

DO 3 1=1, NPAR 
GR AC P ( I )=0.0 
DO ) J=l»N POINT 

) GRACPl I ) =GR ADP ( I H2.0*Z(J, T )*F(J) 

WRITE (6,5007) 

5C0 I FORMAT {//• GRADPHl(l) GRADPHI ( 2 ) GRA 

l DPH I ( 3 ) GRADPHI ( 4 ) GRADPH 1(5)*/) 

WRITE (6,5008) (oRADP( I ), 1*1, NPAR) 

5008 FORMAT (10X,5F18.B> 
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DO 4 1*1 , NPAR 
BETAU )*0. 

DU 4 J* l , NPAR 

BETA* F )*BETA( I >-( L •/?•>*( S ( I , J ) *GRADP ( J ) ) 

PH 1=0. 

DO 5 1*1 , NPO I NT 
5 PH I * PH I »F( I )*F ( I ) 

WRITE (6,5009) 

5C09 FORMAT (//* C(l) C(2) C<3) C(4) C(5) 

1 BETAU) BE T A ( 2 ) BET At 3 ) BETAU) BETAt 5 ) PHI'/) 

WRITE (6,5010) (Cll) , I* l, NPAR) , (BETAt l ) , I * l , NP AR ) , PH I 
5010 FORMAT (1 0 ( F 10 . 5 , l X ) , E9 . 2 > 

WRITE (6,5011) 

5C11 FORMAT (/• EXIT REGRES SUBROUTINE.') 

RETURN 

FND 


SUBROUTINE GRASER ( PH I , C 1 , N , DEI TAC • K X , NSHSO.NPO I NT , NPAR ) 

DIMENSION Cl ( IO),C( 10), CELT AC ( 10) ,DELCMN( 10) 

WRITE (6,5001) 

5C0 l FORMAT ( // * ENTER GRASER SUBROUTINE.') 

00 1 1=1, NPAR 

1 C ( I) =C1 f I ) ♦ DE L T AC t I ) 

CALL SUMSQR ( C , PH I ? , K X , NSMSQ * NPO I N T , NP AR ) 

IF ( K X ) 3,3,? 

2 WRITE (6,5002) 

5002 FORMAT (/' INITIAL CELTAC RESULTS IN AN UNSTABLE SOLUTION. '/' EXIT 
1 GRASER SUBROUTINE.'//) 
u 0 TC 39 

3 DO 4 1 = 1 , NPAR 
DELTACt I ) = DE L T AC ( I )/2. 

4 Ct I >=C1( I ) ♦CELT AC ( I ) 

CALL SUMSQK ( C , PH I 1 , K X , NSMSQ , NPO I N T , NP AR ) 

IF ( K X ) 6,6,5 

5 WRITE (6,5003) 

SCO 3 FORMAT t/' THE STABLE REGION IN PARAMETER SPACE IS NOT CONVEX. •/• 
IEXIT GRASER SUBRUUT I NE . • / / ) 

GO TC 39 

6 IMPHI1-PHI2) 9,9,7 

7 IFIPM2-PHI) 16,8,8 

8 PH I 2* PH I 1 
N = N^ l 

GO TC 3 

9 UCj 10 1 = 1, NPAR 

DE L T AC ( I ) = DE L T AC ( l )/?. 

10 C( 1 )=CU I MDFLTAtl 1 ) 

N*N* 1 

CALL SUMSQR (C,PHlO,KX, NSMSQ, NPOINT ,NPAR ) 

IF(KX) 11,11,5 

11 I F ( PH I 1 — PH I 0 ) 13,13,12 

12 PH 1 2 * PH I 1 
PHI l = PH I 0 
GU TC 9 

13 IF(Phll-PHl) 14, 14, 12 

14 DU 15 1*1, NPAR 

15 DELTACI 1 )=2.*DELTAC( I ) 

GO TC 32 

16 IFIN) 17,17,19 

17 WRITE (6,5004) 

5004 FORMAT (/• THE STEP FOR N*0 IS LOCALLY M I N I M I l I NG . ' / ' EXIT GRASER 
L SUBRGUT I NE • • // ) 

DO 18 1=1, NPAR 
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C l C n»CHl)+ 2 .*DfcLTAC( I ) 
phi*phi3 

GO TG 39 

19 DO 20 1*1, NP AR 

DEL T AC ( I ) =4 • *DE L I AC ( l ) 

20 C ( I ) =C l ( I ) +DELTAC I I ) 

21 N = N-l 

PH I 0 *PH 1 l 
PH I 1 * PH 1 2 

CALL SUMSQR (C,PHI2,KX,NSMSQ,NP0INT,NPAR) 

IF(KX) 23,23,22 

22 WRITE (6,5005) 

5005 FORMAT (/• DECREASING N HAS CAUSED THE SEARCH TO ENTER AN UNSTABLE 
1 REGION. V EXIT GRASER SUBROUTINE.*//) 

GO TC 39 

23 I F ( PH 1 l — PH 12) 30,24,24 

24 IF(NI 25,25,28 

25 PH I = PH I 2 
WRITE 16,5006) 

5006 FORMAT I/' CECRtASING N HAS RtDUCFD N TO ZERO.*/* EXIT GRASER SUBR 
IUUTINF.*//) 

26 00 2 t I = 1 , NP AR 

21 ClIIKKI )*CELTAC( 1 ) 

GD TC 39 

28 DU 29 1 = 1 , NP AR 
DELTACl I )=2.*DFLTAC( I ) 

29 C( 1 )=CL ( I ) ♦DEL T AC ( l ) 

GU TC 21 

30 DO 31 1=1, NP AH 

31 l)E L T AC ( I ) =DE L T AC ( I )/2. 

N = N «■ 1 

32 wUACKN= I 3./4. )* (PH 12-5.* PH I 1 +4 . *PH I 0 ) / ( PH I 2- 3 • *PH | H2. *PHI0) 

DC 33 1=1, N PAR 

DELCFNI I )=CUACMN*DELTAC( I ) 

33 C ( I ) =C l ( I ) ♦ CELCMN ( I ) 

CALL SUPSCR (C,PH IMIN,KX,NSMSQ*NPGINT , NP AR ) 

1 F ( K X ) 34,34,5 

34 IF(PHIPIN-PHU) 35,37,37 

35 00 36 1=1, NP AR 

36 Cl ( I )=C1 ( 1 ) ♦ CELCMN ( I ) 

PHI =PHIMIN 

WRITE (6,5007) PHI 

5C07 FORMAT (/• THE CUACRATIC FIT FORMULA WAS USED TO COMPUTE DELTAC. T 
I HE MINIMUM VALUE FOUND FOR PHI WAS PHI * *,615.8,/* EXIT GRASER SU 
2HR0UT INE. '// ) 

GO TC 39 
3/ PH I = PH I 1 

DO 38 I = 1 , NP AR 

38 Cl I I I =C l (I ) ♦DELTAC( I ) 

WRITE (6,5008) PHI 

5008 FORMAT f/* THE BINARY MINIMUM IS LOWER THAN THE QUADRATIC MINIMUM. 

1 THE MINIMUM VALUE FOUND FOR PHI WAS PHI * *,615.8,/' EXIT GRASER 
2SUBR0UT1NE.'//) 

39 CONTINUE 
RETURN 
END 


SUBROUTINE GKPREX (C,DELTAC,MINPAR,MAXPAR,GRADPfKEXI T , EBDR Y , NP AR ) 
DIMENSION C( 10) , DELTACl 10),GRADP{ 10),C0NBND( 10) 

REAL ko,kstep,minpar< 10).MAXPAR< 10) 

INTEGER CONBND , C I NDE X 
WRITE (6,5001) 
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5001 F ORMAT(/* ENTER GRPREX SUBROUTINE.*) 

C 1 NDEX*0 

K0*0 • 

KEXIT=0 

00 5 I * 1 • NPAR 
CONBNDI I ) =0 

1 F ( C ( I >-M l NP AR ( I) — EBORY) 1*1,3 

1 IF! GRAOP ( I ) ) 5*2*? 

2 OR AUP ( I ) =0. 

C I NOEX=C I NUE X 4 l 
CONflND ( I ) = 1 
GO T(J 5 

3 IF (Cl 1 )-MAXPAR ( I ) ♦EBORY ) 5*4*4 

4 IFIGRAOP(D) 2,2,5 

5 CUNTINUE 

I F I C I NDF X— NPAR ) 8,6,6 

6 WRITE 16*5002) 

5002 FORMAT I/' I HE SEARCH PROCEDURE HAS ATTAINED A CONSTRAINED MINIMUM 
l.*/* EXIT GRPREX SUBROUTINE.*//) 

KE X I T * 1 

00 7 1 = 1 * NPAR 

7 UELTACI I )=0. 

GO IL 17 

0 00 15 1=1, NPAR 

IF (CLNBNDI I ) ) 9,9*15 

9 IF(CRADPIt)) 10,15,11 

10 KSTEP=(C( I ) -MAX PAR I I ) J/GRADPI I ) 

go rc: 12 

11 K S T E P= I C ( 1 ) -MI NPAR ( I ) )/GRADP( I 1 

12 IF (KOI 14, 14, 13 

1 3 IF (KSTFP-KO) 14,15, 15 

14 K0=KSTFP 

15 CUM INUE 

DO 16 1=1, NPAR 

16 OELTACI I ) = — K0*GRADP ( I ) 

WRITF <6,5003) 

5C03 FORMA! <//• DECT AC ( 1 ) DECT AC < 2 ) OE L 

II ACC 3) CELT AC ( 4 ) DEL T AC ( 5 ) • / ) 

WRITE (6,5004) < OELTACI I) , 1= 1, NPAR > 

5004 FORMAT < 10X, 5E10.B) 

WRITE <6,5005) 

5005 FORMAT </• GRADIENT PROJECTION AND EXTRAPOLATION TO A BUUNOARY HAS 
l BEEN ACCGMPL I SHED. • / • EXIT GRPREX SUBROUTINE.*) 

17 CUNTINUE 
RETURN 
END 


SUBROUTINE RANSER < C l , M I NP AR , M A XP AR , MR AND * I Y , NPO I N T , NP AR ) 
DIMENSION Cl < 10), Cl 10), C< 10) 

REAL MINPAR( 10) ,MAXPAR< 10) 

WRITE (6*5001) 

5001 FORMAT ( • l ENT EK RANSER SUBROUTINE.*/) 

NSMSC=0 

UU 1 1=1, NPAR 

1 0(1) =MAXPAR <I)-M INPARU ) 

I X = ! Y 

2 1 = 1 

3 CONTINUE 

CALL RANUUI IX, I Y , YFL ) 

I X = I Y 

Cl (I )=MINPAR< 1)40(1 )*YFL 
T = I 4 I 
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IF <I-NP A R, 3,3,4 

4 CALL SUMSQR ( C 1 , PH I , KX , NSMSQ , NPO I N T , NPAR ) 

IF (KX) 5,5,2 

5 WRITE (6,5002) NSMSC 

5C02 FORMAT (' RANDOM SEARCHING HAS ESTABLISHED A STARTING VALUE FOR PH 
21 . NSMSQ =M3/> 

DC II I s 2 , MR AND 
J* I 

6 CONTINUE 

CALL RANDUt IX, I Y , Y F L ) 

I X = I Y 

C( J)=MINPAK( JI*C<J)*YFL 
J = J* i 

IF ( J-NP AK ) 6,6,7 

7 CALL SUMSOK ( C , PH I R AN , K X , NSMSO* NPO I N T , NP AR ) 

IF ( KX ) H, H , 1 1 

B IF (PHI- PFt IRAN) 11,11,9 

9 PH I = PH I RAN 

UC 10 J=l,NPAR 

10 CMJKIJ) 

11 CONTINUE 

HRITF (6,5003) PHI 

5003 FORMAT I/' THE SMALLEST VALUE FOUNO BY THE RANSER SUBROUTINE FOR T 
lHfc SUM-SQUARED ERROR IS PHI a *E15.8//) 

WRIIF ( 6 , 600 A ) 

5004 FORMAT (• MINIMIZING PARAMETER 

1 VALUES*// 1 C( 1) C ( 2 ) C ( 3 ) 

2 C ( 4 ) C ( 5 ) * / ) 

rfRIlE (6,5005) (Cl( I ) » I - 1 , NP AA ) 

5C0 5 FORMAT (5X,F10*5,4( 10X,F10*5)) 
hRIIE (6,5006) 

5006 FURMAT (/• tXII RANSER SU0ROU T I N6 • f ) 

return 

END 
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